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Abstract 

Synchronization processes in populations of locally interacting elements are in the focus of intense research 
in physical, biological, chemical, technological and social systems. The many efforts devoted to understand 
synchronization phenomena in natural systems take now advantage of the recent theory of complex net- 
works. In this review, we report the advances in the comprehension of synchronization phenomena when 
oscillating elements are constrained to interact in a complex network topology. We also overview the new 
emergent features coming out from the interplay between the structure and the function of the underly- 
ing pattern of connections. Extensive numerical work as well as analytical approaches to the problem are 
presented. Finally, we review several applications of synchronization in complex networks to different dis- 
ciplines: biological systems and neuroscience, engineering and computer science, and economy and social 
sciences. 
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1. Introduction 



Synchronization, as an emerging phenomenon of a population of dynamically interacting units, has 
fascinated humans from ancestral times. Synchronization processes are ubiquitous in nature and play a very 
important role in many different contexts as biology, ecology, climatology, sociology, technology, or even in 
arts P, 0| • It is known that synchrony is rooted in human life from the metabolic processes in our cells to 
the highest cognitive tasks we perform as a group of individuals. For example, the effect of synchrony has 
been described in experiments of people communicating, or working together with a background of shared, 
non-directive conversation, song or rhythm, or of groups of children interacting to an unconscious beat. In 
all cases the purpose of the common wave length or rhythm is to strengthen the group bond. The lack of 
such synchrony can index unconscious tension, when goals cannot be identified nor worked towards because 
the members are "out of sync" 

Among the efforts for the scientific description of synchronization phenomena, there are several capital 
works that represented a breakthrough in our understanding of these phenomena. In 1665, the mathemati- 
cian and physicist, inventor of the pendulum clock, C. Huygens, discovered an odd "kind of sympathy" in 
two pendulum clocks suspended side by side of each other. The pendulum clocks swung with exactly the 
same frequency and 180 degrees out of phase; when the pendula were disturbed, the antiphase state was 
restored within half an hour and persisted indefinitely. Huygens deduced that the crucial interaction for this 
effect came from "imperceptible movements" of the common frame supporting the two clocks. From that 
time on, the phenomenon got into the focus of scientists. Synchronization involves, at least, two elements 
in interaction, and the behavior of a few interacting oscillators has been intensively studied in the physics 
and mathematics literature. However, the phenomenon of synchronization of large populations is a different 
challenge and requires different hypothesi s to be so lved. We will focus our attention on this last challenge. 

In the obituary of Arthur T. Winfree, IStrogatz ^ summarizes what can be consider ed the b eginning of 



the modern quest to explain the synchronization of a population of interacting units: jWienen [5| posed a 
problem in his book Cybernetics: How is it that thousands of neurons or fireflies or crickets can suddenly 
fall into step with one another, all firing or flashing or chirping at the same time, without any leader 
or signal from the environment? Wien er did no t make significant mathematical progress on it, nor did 
anyone else until Winfree came along" . Winfreel ^ studied the nonlinear dynamics of a large population 
of weakly coupled limit-cycle oscillators with intrinsic frequencies that were distributed about some mean 
value, according to some prescribed probability distribution. The milestone here was to consider biological 
oscillators as phase oscillators^ neglecting the amplitude. Working within the framework of a mean field 
model, Winfree discovered that such a population of non-identical oscillators can exhibit a remarkable 
cooperative phenomenon. When the variance of the frequencies distribution is large, the oscillators run 
incoherently, each one near its natural frequency. This behavior remains when reducing the variance until a 
certain threshold is crossed. However, below the threshold the oscillators begin to synchronize spontaneously 
(see 0). Note that the original Winfree model was not solved analytically until recently Q. 

Although Winfree's approach proved to be successful in describing the emergence of spontaneous order 
in the system, it was based on the premise that every oscillator feels the same pattern of interactions. 
However, this all-to-all connectivity between elements of a large population is difficult to conceive in real 
world. When the number of elements is large enough, this pattern is incompatible with physical constraints 
as for example minimization of energy (or costs), and in general with the rare observation of long range 
interactions in systems formed by macroscopic elements. The particular local connectivity structure of the 
elements was missing (in fact, discarded) in these and subsequent approaches. 

In 1998, Watts and Strogatz presented a simple model of network structure, originally intended precisely 
to introduce the connectivity substrate in the problem of synchronization of cricket chirps, which show a high 
degree of coordination over long distances as though the insects were "invisibly" connected. Remarkably, 
this work did not end m a new contribution to synchronization theory but as the seed for the modern theory 
of complex networks [sj. Starting with a regular lattice, they showed that adding a small number of random 
links reduces the distance between nodes drastically, see Fig. [TJ This feature, known as small-world (SW) 
effect, had been first reported in an experiment conducted by S. Milgram [10] examining the average path 
length for social networks of people in the United States. Nowadays, the phenomenon has been detected 
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Figure 1: Smal l-world network cons truction from a regular lattice by rewiring links with a certain probability (randomness), 
as proposed by I Watts and Strogatj [^. 

in many other natural and artificial networks. The inherent complexity of the new model, from now on 
referred to as the Watts-Strogatz (WS) model, was in its mixed nature in between regular lattices and 
random graphs. Very soon, it turned out that the nature of many interaction patterns observed in scenarios 
as diverse as the Internet, the World-Wide Web, scientific collaboration networks, biological networks, was 
even more "complex" than the WS model. Most of them showed a heavy tailed distribution of connectivities 
with no characteristic scale. These networks have been since then called scale-free (SF) networks and the 
most connected nodes are called hubs. This novel structural complexity provoked an explosion of works, 
mainly from the physicists community, since a completely new set of measures, models, and techniques, was 
needed to deal with these topological structures. 

During one decade we have witnessed the evolution of the field of complex networks, mainly from a static 
point of view, although some attempts to characterize the dynamical properties of complex networks have 
also been made. One of these dynamical implications, addressed since the very beginning of the subject, 
is the emergent phenomena of synchronization of a population of units with an oscillating behavior. The 
analysis of synchronization processes has benefited from the advance in the understanding of the topology 
of complex networks, but it has also contributed to the understanding of general emergent properties of 
networked systems. The main goal of this review is precisely to revise the research undertaken so far in 
order to understand how synchronization phenomena are affected by the topological substrate of interactions, 
in particular when this substrate is a complex network. 

The review is organized as follows. We first introduce the basic mathematical descriptors of complex 
networks that will be used henceforth. Next, we focus on the synchronization of populations of oscillators. 
Section IV is devoted to the analysis of the conditions for the stability of the fully synchronized state using 
the Master Stability Function (MSF) formalism. Applications in different fields of science are presented 
afterwards and some perspectives provided. Finally, the last section rounds off the review by giving our 
conclusions. 

2. Complex networks in a nutshell 

There exist excellent reviews devoted to the structural characterization and evolution of complex networks 
11, 12, [l^ 14, [l^, [l^. Here we summarize the main features and standard measures used in complex 



networks. The goal is to provide the reader a brief overview of the subject as well as to introduce some 
notation that will be used throughout the review. 

The mathematical abstraction of a complex network is a graph Q comprising a set of N nodes (or 
vertices) connected by a set of M links (or edges), being ki the degree (number of links) of node i. This 
graph is represented by the adjacency matrix A, with entries = 1 if a directed link from j to i exists, 
and otherwise. In the more general case of a weighted network, the graph is characterized by a matrix 
W, with entries Wij, representing the strength (or weight) of the link from j to i. The investigation of the 
statistical properties of many natural and man-made complex networks revealed that, although representing 
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very different systems, some categorization of them is possible. The most representative of these properties 
refers to the degree distribution P{k), that indicates the probabihty of a node to have a degree k. This 
fingerprint of complex networks has been taken for a long time as its most differentiating factor. However, 
several other measures help to precise the categorization. Examples are the average shortest path length 
£ — (dij), where dij is the length of the shortest path between node i and node j, and the clustering 
coefficient C that accounts for the fraction of actual triangles (three vertices forming a loop) over possible 
triangles in the graph. 

The first classification of complex networks is related to the degree distribution P{k). The differentiation 
between homogeneous and heterogeneous networks in degree is in general associated to the tail of the 
distribution. If it decays exponentially fast with the degree we refer to as homogeneous networks, the most 



representative example being the Erdos-Renyi (ER) random graph 17|. On the contrary, when the tail is 
heavy one can say that the network is heterogeneous. In particular, SF networks are the class of networks 
whose distribution is a power-law, P{k) ~ fc""*", the Barabasi- Albert (BA) model [3| being the paradigmatic 
model of this type of graph. This network is grown by a mechanism in which all incoming nodes are linked 
preferentially to the existing nodes. Note that the limiting case of lattices, or regular networks, corresponds 
to a situation where all nodes have the same degree. 

This categorization can be enriched by the behavior of £. For a lattice of dimension d containing N 
vertices, obviously, £ ^ N^/'^. For a random network, a rough estimate for I is also possible. If the average 
number of nearest neighbors of a vertex is fc, then about k^ vertices of the network are at a distance £ from 
the vertex or closer. Hence, N ^ k^ and then £ ~ \n{N)/ In(fc) , i.e., the average shortest-path length value 
is small even for very large networks. This smallness is usually referred to as the SW property. Associated 
to distances, there exist many measures that provide information about "centrality" of nodes. For instance, 
one can say that a node is central in terms of the relative distance to the rest of the network. One of the 
most frequently used centrality measures in the physics literature is the betweenness (load in some papers), 
that accounts for the number of shortest paths between any pair of nodes in the network that go through a 
given node or link. 

The clustering coefficient C is also a discriminating property between different types of networks. It is 
usually calculated as follows: 



1 " 1 ^ n- 



where is the number of connections between nearest neighbors of node i, and ki is its degree. A large 
clustering coefficient implies many transitive connections and consequently redundant paths in the network, 
while a low C implies the opposite. 

Finally, it is worth mentioning that many networks have a community structure, meaning that nodes 
are linked together in densely connected groups between which connections are sparser. Finding the best 
partition of a network into communities is a very hard problem. The most successful solutions, in terms of 
accuracy and computational cost , are those based on the optimization of a magnitude called modularity, 
proposed in [l^l, that precisely allows for the comparison of different partitionings of the network. The 
modularity of a given partition is, up to a multiplicative constant, the number of links falling within groups 
minus its expected number in an equivalent network with links placed at random. Given a network parti- 
tioned into communities, the mathematical definition of modularity is expressed in terms of the adjacency 
matrix and the total number of links M ~ ^ - ki as 

where Ci is the community to which node i is assigned and the Kronecker delta function i5ci,c takes the 
value 1 if nodes i and j are in the same community, and otherwise. The larger the Q the more modular 
the network is. This proper ty p romises to be specially adequate to unveil structure-function relationships 



in complex networks [21 



f)er ty p ron 
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3. Coupled phase oscillator models on complex networks 



The need to understand synchronization, mainly in the context of biological neural networks, promoted 
the first studies of synchronization of coup led oscilla tors considering a network of interactions between them. 
In the late 80's, Strogatz and Mirollol 25j and later iNiebur et aL 2^ studied the collective synchronization 



of phase non-linear oscillators with random intrinsic frequencies under a variety of coupling schemes in 2D 
lattices. Beyond the differences with the actual conception of a complex network, the topologies studied 
in [2^ can be thought of as a first approach to reveal how the complexity of the connectivity affects 
synchronization. The authors used a square lattice as a geometrical reference to construct three different 
connectivity schemes: four nearest neighbors, Gaussian connectivity truncated at 2a, and finally a random 
sparse connectivity. These results showed that random long-range connections lead to a more rapid and 
robust phase locking between oscillators than nearest-neighbor coupling or locally dense connection schemes. 
This observation is at the root of the recent findings about synchronization in complex networks of oscillators. 
In the current section we review the results obtained so far on three different kinds of oscillatory ensembles: 
limit cycle oscillators (Kuramoto), pulse-coupled models, and finally coupled map systems. We reserve 
for Sect, m those works that use the MSF formalism. Many other works whose major contribution is the 
understanding of synchronization phenomena in specific scenarios are discussed in the Applications section. 

3.1. Phase oscillators 
3.1.1. The Kuramoto model 

The pioneering work by IWinfree 6l| spurred the field of collective synchronization and called for math- 



ematical approaches to tackle the problem. One of these approaches, as already stated, considers a system 
made up of a huge population of weakly-coupled, nearly identical, interacting limit-cycle oscillators, where 
each oscillator exerts a phase dependent influence on the others and changes its rhythm according to a 
sensitivity function [13, [2^ . 

Even if these simplifications seem to be very crude, the phenomenology of the problem can be captured. 
Namely, the population of oscillators exhibits the dynamic analog to an equilibrium phase transition. When 
the natural frequencies of the oscillators are too diverse compared to the strength of the coupling, they 
are unable to synchronize and the system behaves incoherently. However, if the coupling is strong enough, 
all oscillators freeze into synchrony. The transition from one regime to the other takes place at a certain 
threshold. At this point some elements lock their relative phase and a cluster of synchronized nodes develops. 
This constitutes the onset of synchronization. Beyond this value, the population of oscillators is split into 
a partially synchronized state made up of oscillators locked in phase and a group of nodes whose natural 
frequencies are too different as to be part of the coherent cluster. Finally, after further increasing the 
coupling, more and more elements get entrained around the mean phase of the collective rhythm generated 
by the whole population and the system settles in the completely synchronized state. 

Kuramoto [29, 30] worked out a mathematically tractable model to describe this phenomenology. He 
recognized that the most suitable case for analytical treatment should be the mean field approach. He 
proposed an all-to-all purely sinusoidal coupling, and then the governing equations for each of the oscillators 
in the system are: 

K ^ 

where the factor is incorporated to ensure a good behavior of the model in the thermodynamic limit, 
N —> 00, LUi stands for the natural frequency of oscillator i, and K is the coupling constant. The frequencies 
LUi are distributed according to some function g{uj), that is usually assumed to be unimodal and symmetric 
about its mean frequency fl. Admittedly, due to the rotational symmetry in the model, we can use a rotating 
frame and redefine uji uji + Q for all i and set f2 = 0, so that the uji^s denote deviations from the mean 
frequency. 
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The collective dynamics of the whole population is measured by the macroscopic complex order param- 
eter, 

1 ^ 

Kt)e'^'*^ = ^Ee''^^*^ (4) 

where the modulus < r(t) < 1 measures the phase coherence of the population and 4>{t) is the average 
phase. The values r ~ 1 and r ~ (where ~ stands for fluctuations of size 0(7V~i/2)) describe the limits 
in which all oscillators are either phase locked or move incoherently, respectively. Multiplying both parts of 
Eq. ([4]) by e~'^* and equating imaginary parts gives 



r sm( 



1 ^ 

J2smi9,-d,), (5) 



N 



yielding 

0, = CO, + Kr sin {(j)-e,) {i = l,...,N). (6) 

Equation ^ states that each oscillator interacts with all the others only through the mean field quantities r 
and (f). The first quantity provides a positive feedback loop to the system's collective rhythm: as r increases 
because the population becomes more coherent, the coupling between the oscillators is further strengthened 
and more of them can be recruited to take part in the coherent pack. Moreover, Eq. © allows to calculate 
the critical coupling Kc and to characterize the order parameter limt^oo ft{K) — r{K). Looking for steady 
solutions, one assumes that r{t) and (l>(t) are constant. Next, without loss of generality, we can set = 0, 
which leads to the equations of motion [2^, [s^l 

e, = uj,-Kr sin (Oi) {i^l,...,N). (7) 

The solutions of Eq. ([7]) reveal two different types of long-term behavior when the coupling is larger than 
the critical value, Kc- On the one hand, a group of oscillators for which jcjil < Kr are phase-locked at 
frequency f2 in the original frame according to the equation uJi = Kr sin (9i). On the other hand, the 
rest of the oscillators for which \uji\ > Kr holds, are drifting around the circle, sometimes accelerating 
and sometimes rotating at lower frequencies. Demanding some conditions for the stationary distribution of 
drifting oscillators with frequency uJi and phases 9i [13], a self-consistent equation for r can be derived as 

r = Kr r (cos^ 9) g{uj)dd, 



where u = Kr sin (9). This equation admits a non-trivial solution, 

beyond which r > 0. Equation ([5]) is the Kuramoto mean field expression for the critical coupling at the 
onset of synchronization. Moreover, near the onset, the order parameter, r, obeys the usual square-root 
scaling law for mean field models, namely, 

r^{K~ Kef (9) 

with P = 1/2. Numerical simulations of the model verified these results. The Kuramoto model (KM, from 
now on) approach to synchronization was a breakthrough for the understanding of synchronization in large 
populations of oscillators. 

Even in the simplest case of a mean field interaction, there are still unsolved problems that have resisted 
any analytical attempt. This is the case, e.g., for finite populations of oscillators and some questions 
regarding global stability results In what follows, we focus on another aspect of the model's assumptions, 
namely that of the connection topology of real systems [3, [H, which usually do not show the all-to-all 



pattern of interconnections underneath the mean field approach. 

7 



3.1.2. Kuramoto model on complex networks 

To deal with the KM on complex topologies, it is necessary to reformulate Eq. ^ to include the 
connectivity 



aijOij sin{6j ~ 0i) (« = 1, TV) , 



(10) 



i 



where (Jy is the coupling strength between pairs of connected oscillators and a^- are the elements of the 
connectivity matrix. The original Kuramoto model is recovered by letting a^- = l,Vi 7^ j (all-to-all) and 



The first problem when defining the KM in complex networks is how to state the interaction dynamics 
properly. In contrast with the mean field model, there are several ways to define how the connection topology 
enters in the governing equations of the dynamics. A good theory for Kuramoto oscillators in complex 
networks should be phenomenologically relevant and provide formulas amenable to rigorous mathematical 
treatment. Therefore, such a theory should at least preserve the essential fact of treating the heterogeneity 
of the network independently of the interaction dynamics, and at the same time, should remain calculable 
in the thermodynamic limit. 

For the original model, Eq. ([3]), the coupling term on the right hand side of Eq. pil| is an intensive 
magnitude because the dependence on the size of the system cancels out. This independence on the number of 
oscillators TV is achieved by choosing cr^ = K/N . This prescription turns out to be essential for the analysis 
of the system in the thermodynamic limit N ^ 00 ixi the all-to-all case. However, choosing aij = K/N 
for the governing equations of the KM in a complex network makes them to become dependent on N . 
Therefore, in the thermodynamic limit, the coupling term tends to zero except for those nodes with a degree 
that scales with N . Note that the existence of such nodes is only possible in networks with power-law degree 
distributions ITsj . but this happens with a very small probability as k with 7 > 2. In these cases, 
mean field solutions independent of N are recovered, with slight differences in the onset of synchronization 
of all-to-all and SF networks [31] . 

A second prescription consists in taking — K/ki (where ki is the degree of node i) so that is a 
weighted interaction factor that also makes the right hand side of Eq. (jTU]) intensive. This form has been 
used to solve the paradox of heterogeneity [s^ that states that the heterogeneity in the degree distribution, 
which often reduces the average distance between nodes, may suppress synchronization in networks of 
oscillators coupled symmetrically with uniform coupling strength. This result refers to the stability of the 
fully synchronized state, but not to the dependence of the order parameter on the coupling strength (where 
partially synchronized and unsynchronized states exist). Besides, the inclusion of weights in the interaction 
strongly affects the original KM dynamics in complex networks because it can impose a dynamic homogeneity 
that masks the real topological heterogeneity of the network. 

The prescription aij = /C/const, which may seem more appropriate, also causes some conceptual prob- 
lems because the sum in the right hand side of Eq. (|10|) could eventually diverge in the thermodynamic 
limit. The constant in the denominator could in principle be any quantity related to the topology, such as 
the average connectivity of the graph, (A:), or the maximum degree A:max- Us physical meaning is a re-scaling 
of the temporal scales involved in the dynamics. However, except for the case of aij = i^/fcmax, the other 
possible settings do not avoid the problems when — > 00. On the other hand, for a proper comparison of 
the results obtained for different complex topologies (e.g. SF or uniformly random), the global and local 
measures of coherence should be represented according to their respective time scales. Therefore, given two 
complex networks A and B with /cmax — kA and fcmax = fcs respectively, it follows that to make meaningful 
comparisons between observables, the equations of motion Eq. (fTU)) should refer to the same time scales, 
i.e., aij = KA/kA = Ks/ks = a. With this formulation in mind, Eq. (|10|) reduces to 



independently of the specific topology of the network. This allows us to study the dynamics of Eq. (fTT|) on 
different topologies, compare the results, and properly inspect the interplay between topology and dynamics 
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K/N, \fi,j. 





(11) 



in what concerns synchronization. 

As we shall see, there are also several ways to define the order parameter that characterizes the global 
dynamics of the system, some of which were introduced to allow for analytical treatments at the onset of 
synchronization. We advance, however, that the same order parameter, Eq. (j4]), is often used to describe 
the coherence of the synchronized state. 



3.1.3. Onset of synchronization in complex networks 

Studies on synchronization in complex top olog ies where each node is considered to be a Kuramoto 
oscillator, were first reported for WS networks [33, [31 and BA graphs 'si']. These works are mainly 
numerical explorations of the onset of synchronization, their main goal being the characterization of the 



critical coupling beyond which groups of nodes beating coherently first appear. In 3J|, the authors considered 
oscillators with intrinsic frequencies distributed according to a Gaussian distribution with unit variance 
arranged in a WS network with varying rewiring probability, p, and explored how the order parameter, 
Eq. ([1]), changes upon addition of long-range links. Moreover, they assumed a normalized coupling strength 
Uij — K/{k) , where (fc) is the average degree of the graph. Numerical integration of the equations of motion 
(|10|1 under variation of p shows that collective synchronization emerges even for very small values of the 
rewiring probability. 

The results confirm that networks obtained from a regular ring by just rewiring a tiny fraction of links 
ip ^ 0) can be synchronized with a finite K. Moreover, in contrast with the arguments provided in [S^l, 
we notice that their results had been obtained for a fixed average degree and thus the Kuramoto's critical 
coupling can not be recovered by simply taking p —^ 1, which produces a random ER graph with a fixed 
minimum connectivity. This limit is recovered by letting (fc) increase. Actually, numerical simulations of 
the same model in [33| showed that the Kuramoto limit is approached when the average connectivity grows. 

In [ssl the same problem in BA networks it is considered. The natural frequencies and the initial values 
of 0i were randomly drawn from a uniform distribution in the interval (—1/2, 1/2) and (— tt, tt), respectively. 
The global dynamics of the system, Eq. (jlip , turns out to be qualitatively the same as for the original KM 
as shown in Fig. [21 where the dependence of the order parameter Eq. ^ with a is shown for several system 
sizes. 

The existence of a critical point for the KM on SF networks came as a surprise. Admittedly, this is 
one of the few cases in which a dynamical process shows a critical behavior when the substrate is described 
by a power-law connectivity distribution with an exponent 7 < 3 14, 37 1. In principle it could be a 
finite size effect, but it turned out from numerical simulations that this was not the case. To determine the 
exact value of CTc, one can make use of standard finite-size scaling analysis. At least two complementary 
strategies have been reported. The first one allows bounding the critical point and is computationally more 
expensive. Consider a network of size N, for which no synchronization is attained below CTc, where r{t) 
decays to a small residual value of size 0(1/ ^/N). Then, the critical point may be found by examining the 
iV-dependence of r{a,N). In the sub-critical regime (cr < CTc), the stationary value of r falls off as N^^^"^, 
while for a > the order parameter reaches a stationary value as — > 00 (though still with 0(l/\/]V) 
fluctuations). Therefore, plots of r versus N allow us to locate the critical point CTc. Alternatively, a more 
accurate approach can be adopted. Assume the scaling form for the order parameter [s^ : 



r = N-^!{N^{a~G,)), (12) 

where /(x) is a universal scaling function bounded as x ^ ±cxd and a and v are two critical exponents to 
be determined. Since at ct = (Tc, the value of the function / is independent of TV, the estimation of CTc can 
be done by plotting A^"r as a function of a for various sizes N and then finding the value of a that gives a 
well-defined crossing point, the critical coupling CTc. As a by-product, the method also allows us to calculate 
the two scaling exponents a and z^, the latter can be obtained from the equality 

\i\{(drlda)\„^\ (i^ - a) IniV + const, (13) 

once a is computed. 
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Figure 2: Order parameter r (Eq. l(4j) as a function of a for several BA networks of different sizes. Finite size scaling analysis 
shows that the onset of synchronization takes place at a critical value = 0.05(1). The inset is a zoom around CTc. From [35ll . 



Following these scaling procedures, it was estimated a value for the critical coupling strength CTc = 0.05(1) 
[35L I39I . l40| . Moreover, r ~ (cr — CTc)^ when approaching the critical point from above with /? = 0.46(2) 
indicating that the square-root behavior typical of the mean field version of the model (/3 = 1/2) seems to 
hold as well for BA networks. 

Before turning our attention to some theoretical attempts to tackle the onset of synchronization, it is 
worth to briefly summarize other numerical results that have explored how the critical coupling depends 
on other topological features of the underlying SF graph. Recent results have shed light on the influence 
of the topology of the local interactions on the route to and the onset of synchronization. In particular, 
the authors in 41, 4^, 43 1 explored the Kuramoto dynamics on networks in which the degree distribution 
is kept fixed, while the clustering coefficient (C) and the average path length (€) of the graph change. The 
results suggest that the onset of synchronization is mainly determined by C, namely, networks with a high 
clustering coefficient promote synchronization at lower values of the coupling strength. On the other hand, 
when the coupling is increased beyond the critical point, the effect of I dominates over C and the phase 
diagram is smoothed out (a sort of stretching), delaying the appearance of the fully synchronized state as 
the average shortest path length increases. 

In a series of recent papers 



311 . |44| . I45l |46| . 1471 . 148j , the onset of synchronization in large networks of 



coupled oscillators has been analyzed from a theoretical point of view under certain (sometimes strong) 
assumptions. Despite these efforts no exact analytical results for the KM on general complex networks are 
available up to date. Moreover, the analytical approaches predict that for uncorrelated SF networks with 
an exponent 7 < 3, the critical coupling vanishes as iV ^ cxd, in contrast to numerical simulations on BA 
model networks. It appears that the strong heterogeneity of real networks and the finite average connectivity 
strongly hampers analytical solutions of the model. 

iji. Defining 



Following [sij, consider the system in Eq. (llip . with a symmetric0 adjacency matrix 



-"^The reader can find the extension of the forthcoming formalism to directed networks in [4J 
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a local order parameter as 

TV 

r,e'^' =^ay(e"'^)t , (14) 

where (• • ■)t stands for a time average, a new global order parameter to measure the macroscopic coherence 
is readily introduced as 

T: 

(15) 



Now, rewriting Eq. pip as a function of r^, yields, 

9i = Lu, - an sm{6i - 4>i) - ahi{t) . (16) 

In Eq. US]), hi{t) = Im{e~'^* ^^j^ ay((e'^j)t — e'^^)} depends on time and contains time fluctuations. 
Assuming the terms in the previous sum to be statistically independent, hi(t) is expected to be proportional 
to ^/ki above the transition, where Vi ^ ki. Therefore, except very close to the critical point, and assuming 
that the number of connections of each node is large enougljl [ki ^ 1 as to be able to neglect the time 
fluctuations entering /ij, i.e., hi <C r^), the equation describing the dynamics of node i can be reduced to 

di = uji ~ art sin(6li - (f>,) . (17) 

Next, we look for stationary solutions of Eq. pT)) . i.e. sin(6',i — (pi) — LOi/ari. In particular, oscillators whose 
intrinsic frequency satisfies jwij < ari become locked. Then, as in the Kuramoto mean field model, there 
are two contributions (though in this case to the local order parameter), one from locked and the other from 
drifting oscillators such that 

N 

n = ^a.,(e'(^^-*')),= (18) 

|ci;j|<crrj |cjj|>crrj 

To move one step further, some assumptions are needed. Consider a graph such that the average degree of 
nearest neighbors is high (i.e., if the neighbors of node i are well-connected). Then it is reasonable to assume 
that these nodes are not affected by the intrinsic frequency of i. This is equivalent to assume solutions (r.^, (j)i) 
that are, in a statistical sense, independent of the natural frequency uoi. With this assumption, the second 
summand in Eq. psp can be neglected. Taking into account that the distribution g{Lij) is symmetric and 
centered at f2 = 0, after some algebra one is left with [sij 



ri= ^ fly cos((/)j - (/)^)W 1 - f ^ j . (19) 

The critical coupling ctc is given by the solution of Eq. that yields the smallest a. It can be argued that 
it is obtained when cos(0j- — 0i) = 1 in Eq. HH), thus 



2 



\u},\<<yrj * \ Jy 

which is the main equation of the time average approximation (recall that time fluctuations have been 
neglected). Note, however, that to obtain the critical coupling, one has to know the adjacency matrix as 



^This obviously restricts the range of real networks to which the approximation can be applied. 
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well as the particular values of LOi for all i and then solve Eq. ([^0]) numerically for the {ri}. Finally, the 
global order parameter defined in Eq. psp can be computed from r^. 

Even if the underlying graph satisfies the other aforementioned topological constraints, it seems unre- 
alistic to require knowledge of the {cjij's. A further approach, referred to as the frequency distribution 
approximation can be adopted. According to the assumption that ki ^ 1 for all i, or equivalently, that the 
number of connections per node is large (a dense graph), one can also consider that the natural frequencies 
of the neighbors of node i follows the distribution g{'jj). Then, Eq. (PO]) can be rewritten avoiding the 
dependence on the particular realization of {tOi} to yield. 



2 



= cr^^ayTj / g{xarj)\/\ — x^dx , (21) 
j •'-^ 

with X — uj/[arj). This equation allows us to readily determine the order parameter r as a function of 
the network topology (a^), the frequency distribution ((/(u;)) and the control parameter (cr). On the other 
hand, Eq. ([^D|) still does not provide explicit expressions for the order parameter and the critical coupling 
strength. To this end, one introduces a first-order approximation g(xarj) « ^(0) which is valid for small, 
but nonzero, values of r. Namely, when rj 0+ 



\ ' 



where Kc ~ 2/{Trg{0)) is Kuramoto's critical coupling. Moreover, as the smallest value of a corresponds to 
CTc, it follows that the critical coupling is related to both Kc and the largest eigenvalue Amax of the adjacency 
matrix, yielding 

ac = ^ . (22) 

Equation states that in complex networks, synchronization is first attained at a value of the coupling 
strength that inversely depends on g{0) and on the largest eigenvalue Amax of the adjacency matrix. Note 
that this equation also recovers Kuramoto's result when mj = 1, ^ j, since Amax = — 1. It is 
worth stressing that although this method allows us to calculate ctc analytically, it fails to explain why for 
uncorrelated random SF networks with 7 < 3 and in the thermodynamic limit N 00, the critical value 
remains finite. This disagreement comes from the fact that in these SF networks, Amax is proportional 

to the cutoff of the degree distribution, fcmax which in turn scales with the system size. Putting the two 

i 1 

dependencies together, one obtains Amax ~ fcmax ~ A?^ 2(7-1) ^ _> go, thus predicting CTc = in 

the thermodynamic limit, in contrast to finite size scaling analysis for the critical coupling via numerical 
solution of the equations of motion. Note, however, that the difference may be due to the use of distinct 
order parameters. Moreover, even in the case of SF networks with 7 > 3, Amax still diverges when we take 
the thermodynamic limit, so that cTc as well. As we shall see soon, this is not the case when other 
approaches are adopted, at least for 7 > 3. 

It is possible t o g o beyond with the latter approximation and to determine the behavior of r near the 
critical point. In |3l| a perturbative approach to higher orders of Eq. ()2ip is developed, which is valid for 



relatively homogeneous degree distributions (7 > 5)|^ They showed that for (ct/ctc) — 1^0+ 



^The approach holds if the fourth moment of the degree distribution, (fc*) = P(k)k'^dk remains finite when N 00. 
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where ryi = —Trg" {0)Kc/l6 and 



2 

max 



(24) 



where u is the normalized eigenvector of the adjacency matrix corresponding to Amax and (u^) — Uj/N. 
The analytical insights discussed so far can also be reformulated in terms of a mean field approximation 



46l . |48| . 1 3X1 . |47| for complex networks. This approach (valid for large enough (k) ) considers that every 



oscillator is influenced by the local field created in its neighborhood, so that is proportional to the degree 
of the nodes ki, i.e., Vi ki. Assuming this is the case and introducing the order parameter r through 



1 



N 



(25) 



after summing over i and substituting — rki in Eq. (j2ip we obtain (sij 

N N 1 

^^fcj=(T^^fc| J g{x(jrkj)\/l — x'^dx 



(26) 



The above relation, Eq. was independently derived in ^io"], who first studied analytically the problem 
of synchronization in complex networks, though using a different approach. Taking the continuum limit, 
Eq. becomes 

J kP{k)dk = a j k^P{k)dk J g{xark)^l-x^dx , 

which for r — > 0"*" verifies 



(27) 



kP(k)dk 



a / k'^P{k)dk / g{Q)y/l-x'^dx 



o-.9(0)7r , ,2 



k^P{k)dk , 



(28) 



which leads to the condition for the onset of synchronization (r > 0) as 

crg(0)7r 



k^P{k)dk > J kP{k)dk 



that is. 



(fc) 



(k) 



7rg(0)(fc2) '"'- {k^) ■ ^^^^ 

The mean field result, Eq. (|29p. gives as a surprising result that the critical coupling CTc in complex networks 
is nothing else but the one corresponding to the all-to-all topology Kc re-scaled by the ratio between 
the first two moments of the degree distribution, regardless of the many differences between the patterns 
of interconnections. Precisely, it states that the critical coupling strongly depends on the topology of 
the underlying graph. In particular, (Tc ^ when the second moment of the distribution (k^) diverges, 
which is the case for SF networks with 7 < 3. Note, that in contrast with the result in Eq. ([2^ . for 
7 > 3, the coupling strength does not vanish in the thermodynamic limit. On the other hand, if the 
mean degree is kept fixed and the heterogeneity of the graph is increased by decreasing 7, the onset of 
synchronization occurs at smaller values of CTc- Interestingly enough, the dependence gathered in Eq. ([^ 
has the same functional form for the critical points of other dynamical processes such as percolation and 
epidemic spreading processes [ij, |l5|, |37| . While this result is still under numerical scrutiny, it would imply 
that the critical properties of many dynamical processes on complex networks are essentially determined by 
the topology of the graph, no matter whether the dynamics is nonlinear or not. The corroboration of this 
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last claim will be of extreme importance in physics, probably changing many preconceived ideas about the 
nature of dynamical phenomena. 

Within the mean field theory, it is also possible to obtain the behavior of the order parameter r near 
the transition to synchronization. Equation (|29p was also independently derived in [48| starting from the 
differential equation Eq. pT|) . Using the weighted order parameter 

and assuming the same magnitude of the effective field of each pair of coupled oscillators one obtains 

Oi = - -^hfsm{e.i) , (30) 

where we have set = 0. Now, it is considered again that in the stationary state the system divides into 
two groups of oscillators, which are either locked or rotating in a nonuniform manner. Following the same 
procedure employed in all the previous derivations, the only contribution to r comes from the former set of 
oscillators. After some algebra it is shown that the critical coupling CTj. is given by Eq. and that 
near criticality 

(31) 

for 7 > 3, with a critical exponent /3 = ^ if 7 > 5, and /? = when 3 < 7 < 5. For the most common 
cases in real networks of 2 < 7 < 3, the critical coupling tends to zero in the thermodynamic limit so that 
r should be nonzero as soon as a =/= 0. In this case, one gets r ^ cr^/^^"''''. Notably, the latter equation is 
exactly the same found for the absence of critical behavior in the region 2 < 7 < 3 for a model of epidemic 
spreading [iot . 

One recent theoretical study in [50| is worth mentioning here. They have extended the mean field 
approach to the case in which the coupling is asymmetric and depends on the degree. In particular, they 
studied a system of oscillators arranged in a complex topology whose dynamics is given by 



N 

,1-1-1 

j=i 



%). (32) 



T) = 1 corresponds to the symmetric, non-degree dependent, case. Extending the mean field formalism to the 
cases rj 1, they investigated the nature of the synchronization transition as a function of the magnitude 
and sign of the parameter rj. By exploring the whole parameter space {rj, a), they found that for rj ~ and 
SF networks with 2 < 7 < 3, a finite critical coupling <Tc is recovered in sharp contrast to the non-weighted 
coupling case for which we already know that ac = 0. This result seems phenomenologically meaningful, 
since setting 77 = implies that the coupling in Eq. (fT(I|) is o-y = a/ki, which, as discussed before [il], might 
have the effect of partially destroying the heterogeneity inherent to the underlying graph by normalizing all 
the contributions X^jLi '^ij sin{9j — 9i) by ki = X^jLi ^ij- 

3.1. 4-. Path towards synchronization in complex networks 

Up to now, we have discussed both numerically and theoretically the onset of synchronization. In the 
next section, we shall also discuss how the structural properties of the networks influence the stability of 
the fully synchronized state. But, what happens in the region where we are neither close to the onset of 
synchronization nor at complete synchronization? How is the latter state attained when different topologies 
are considered? 

As we have seen, the influence of the topology is not only given by the degree distribution, but also by 
how the oscillators interact locally. To reduce the number of degrees of freedom to a minimum, let us focus 
on the influence of heterogeneity and study the evolution of synchronization for a family of complex networks 
which are comparable in their clustering, average distance and correlations so that the only difference is due 
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Table 1: Topological properties of the networks and critical coupling strengths derived from a finite size scaling analyses, Eq. 
Ijl2|l . Different values of x corresponds to grown networks whose degree of heterogeneity varies smoothly between the two 
limiting cases of ER and SF graphs. From 14011 . 



X 




(P) 


h 




0.0 


(SF) 


115.5 


326.3 


0.051 


0.2 




56.7 


111.6 


0.066 


0.4 




44.9 


47.7 


0.088 


0.6 




41.1 


25.6 


0.103 


0.8 




39.6 


16.8 


0.108 


1.0 


(ER) 


39.0 


14.8 


0.122 



to the degree distribution^ For these networks, the previous theoretical approaches argued that the critical 
couphng (Tc is proportional to {k)/{k^) , so that different topologies should give rise to distinct critical 
points. In particular, in [39l l40j it was studied numerically the path towards synchronization in ER and SF 
networks. Theyalso studied several networks whose degree of heterogeneity can be tuned between the two 
limiting cases [5l|- These authors put forward the question: How do SF networks compare with ER ones 
and what are the roots of the different behaviors observed? 

Numerical simulations 3^ iO] confirm qualitatively the theoretical predictions for the onset of synchro- 
nization, as summarized in Table [TJ In fact, the onset of synchronization first occurs for SF networks. 
As the network substrate becomes more homogeneous, the critical point ctc shifts to larger values and the 
system seems to be less synchronizable. On the other hand, they also showed that the route to complete 
synchronization, r = 1, is sharper for homogeneous networks. No critical exponents for the behavior of r 
near the transition points have been reported yet for the ER network, so that comparison with the mean 
field value (3 = 1/2 for a SF network with 7 = 3 is not possibleU Numerically, a detailed finite size scaling 
analysis in SF and ER topologies shows that the critical coupling strength corresponds in SF networks to 
crf^ = 0.051, and in random ER networks to cr™ = 0.122, a fairly significant numerical difference. 

The mechanisms behind the differences in the emergence of collective behavior for ER and SF topologies 
can be explored numerically by defining a local order parameter that captures and quantifies the way in 
which clusters of locked oscillators emerge. The main difference with respect to r is that one measures the 
degree of synchronization of nodes (r) with respect to the average phase </> and the other (runk) to the degree 
of synchronization between every pair of connected nodes. Thus, riink gives the fraction of all possible links 
that are synchronized in the network as 



''link 



1 ^ 
= > a. 



lim — — 

At^oo At 



(33) 



beiiig tr the time the system needs to settle into the stationary state, and At a large averaging time. In 
3^ 40j the degree of synchronization of pairs of connected oscillators was measured in terms of the symmetric 
matrix 



lim — — 

At^oo At 



t,. 



(34) 



which, once filtered using a threshold T such that the fraction of synchronized pairs equals riink, allows us 
to identify the synchronized links and reconstruct the clusters of synchrony for any value of cr, as illustrated 
in FiglSl From a microscopic analysis, it turns out that for homogeneous topologies, many small clusters of 



*This isolation of individual features of complex networks is essential to understand the interplay between topology and 
dynamics. As we will discuss along the review, many times this aspect has not been properly controlled raising results that 
are confusing, contradictory or even incorrect. 

^The numerical value of (3 contradicts the prediction of the mean-field approach (sec the discussion after Eq. I lijll l) The 
reason of such discrepancy is not clear yet. 
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Figure 3: (color online) Synchronized components for several values of a for the two limiting cases of ER and SF networks. 
The figure clearly illustrates the differences in forming synchronization patterns for both types of networks: in the SF case 
links and nodes are incorporated together to the largest of the synchronized clusters, while for the ER network, what is added 
are links between nodes already belonging to such cluster. From [soll . 



synchronized pairs of oscillators are spread over the graph and merge together to form a giant synchronized 
cluster when the effective coupling is increased. On the contrary, in heterogeneous graphs, a central core 
containing the hubs first comes up driving the evolution of synchronization patterns by absorbing small 
clusters. Moreover, the evolution of runk as a grows explains why the transition is sharper for ER networks: 
nodes are added first to the giant synchronized cluster and later on the links among these nodes that were 
missing in the original clusters of synchrony. In SF graphs, oscillators are added to the largest synchronized 
component together with most of their links, resulting in a much slower growth of runk- Finally, it is also 
computed the probability that a node with degree k belongs to the largest synchronized cluster and reported 
that this probability is an increasing function of k for every a, namely, the more connected a node is, the 
more likely it takes part in the cluster of synchronized links [s^, . It is interesting to mention here that 
a similar dependence is obtained if one analyzes the stability of the synchronized state under perturbations 
of nodes of degree k. In [s^ it was found that the average time (r) a node needs to get back into the fully 
synchronized state is inversely proportional to its degree, i.e., (r) ^ k^^. 

Very recently [H^], the path towards synchronization was also studied looking for the relation between 
the time needed for complete synchronization and the spectral properties of the Laplacian matrix of the 
graph, 

Lij kiSij - a.ij. (35) 

The Laplacian matrix is symmetric with zero row-sum and hence all the eigenvalues are real and non- 
negative. Considering the case of identical Kuramoto oscillators, whose dynamics has only one attractor, 
the fully synchronized state, they found that the synchronization time scales with the inverse of the smallest 
nonzero eigenvalue of the Laplacian matrix. Surprisingly, this relation qualitatively holds for very different 
networks where synchronization is achieved, indicating that this eigenvalue alone might be a relevant topo- 



logical property for synchronization phenomena. The authors in [53| remark the role of this eigenvalue not 



only for synchronization purposes but also for the flow of random walkers on the network. 

3.1.5. Kuramoto model on structured or modular networks 

In this section, we discuss a context in which synchronization has turned out to be a relevant phenomenon 
to explore the relation between dynamical and topological properties of complex networks. Many complex 
networks in nature are modular, i.e. composed of certain subgraphs with differentiated internal and external 
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connectivity that form communities |20l . Il9l . |l5| . This is a hmiting situation in which the local structure 
may greatly affect the dynamics, irrespective of whether or not we deal with homogeneous or heterogeneous 
networks. 

Synchronization processes on modular networks have been recently studied as a mechanism for commu- 
nity detection (H^jHH, 5^ 57 1. The situation in which a set of identical Kuramoto oscillators (i.e., coi 
with random initial conditions evolves after a transient to the synchronized state was addressed In 
Note that in this case full synchronization is always achieved as this state is the only attractor of the dy- 
namics so that the coupling strength sets the time scale to attain full synchronization: the smaller a is, the 
longer the time scale. The authors in [55] guessed that if high densely interconnected motifs synchronize 
more easily than those with sparse connections [s^, then the synchronization of complex networks with 
community structure should behave differently at different time and spatial scales. In synthetic modular 
networks, starting from random initial phases, the highly connected units forming local clusters synchronize 
first and later on, in a sequential process, larger and larger topological structures do the same up to the 
point in which complete synchronization is achieved and the whole population of oscillators beat at the 
same pace. This process occurs at different time scales and the dynamical route towards the global attractor 
reveals the topological structures that represent communities, from the microscale at very early states up 
to the macroscale at the end of the time evolution. 

The authors studied the time evolution of pairs of oscillators defining the local order parameter 

p,,it) = {cos[0,{t)-0j{t)]) , (36) 

averaged over different initial conditions, which measures the correlation between pairs of oscillators. To 
identify the emergence of compact clusters reflecting communities, a binary dynamic connectivity matrix is 
introduced such that 

ff p.,it)<T, (3^) 

for a given threshold T. Changing the threshold T at fixed times reveals the correlations between the dynam- 
ics and the underlying structure, namely, for large enough T, one is left with a set of disconnected clusters 
or communities that are the innermost ones, while for smaller values of T inter-community connections 
show up. In other words, the inner community levels are the first to become synchronized, subsequently 
the second level groups, and finally the whole system shows global synchronization. Note that since the 
function Pij(t) is continuous and monotonic, we can define a new matrix Vxit), that takes into account the 
time evolution for a fixed threshold. The evolution of this matrix unravels the topological structure of the 
underlying network at different time scales. In the top panels of Fig. [3] we plot the number of connected 
components corresponding to the binary connectivity matrix with a fixed threshold as a function of time for 
networks with two hierarchical levels of communities. There we can notice how this procedure shows the 
existence of two clear time scales corresponding to the two topological scales. 

It is also possible to go one step further and show that the evolution of the system to the global attractor 
is intimately linked to the whole spectrum of the Laplacian matrix The bottom panels of Fig. H] 

show the ranked index of the eigenvalues of Lij versus their inverse. As can be seen, both representations 
(top and bottom) are qualitatively equivalent, revealing the topological structure of the networks. The only 
difference is that one comes from a dynamical matrix and the other from the spectrum of the Laplacian, 
that fully characterizes the topology. Thus, synchronization can be used to unveil topological scales when 
the architecture of the network is unknown. 

The relationship between the eigenvalue spectrum of Lij and the dynamical structures of Fig. [H can be 



understood from the linearized dynamics of the Kuramoto model, which reads |55l . |5; 

N 



aYUjOj i = l,...,iV, (38) 



It is worth stressing here that for this purpose the assumption oi Wi = u) can be adopted without loss of generality as it 
makes the analysis easier. Synchronization of non-identical oscillators also reveals the existence of community structures. See 
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Figure 4; (color online) Top: Number of disconnected synchronized components (d.s.c.) as a function of time. Bottom: Rank 
index versus the corresponding eigenvalue of the Laplacian matrix. Each column corresponds to a network with two hierarchical 
levels of communities. The difference lies in the relative weight of the two modular levels. From [55| . 



that is a good approximation after a fast transient starting from random initial phases in the range [0, 27r]. 
The solution of Eq. ([55)) , in terms of the normal modes tpi (t) , is 

JV JV 

e,{t) = ^ B,,^/j,{t) = ^ B,,^,(0)e-'^^^* , (39) 
i=i i=i 

where Bij is the normalized eigenvector matrix and \i the eigenvalues of the Laplacian matrix. Going back 
to the original coordinates, the phase difference between any pair of oscillators is 

JV 

\m ~ < - \^]\ I ^'^■(0)1 ■ (40) 

Assuming a global bound for the initial conditions |0fc(O)| < 0, Vfc and taking into account the normal- 
ization of the eigenvector matrix, \Bij\ < 1, we can sum over the index k to get 

JV 

\<^^{t) - e^{t)\ < A^e^ |By - B.mj\ e-"^^*. (41) 

The sum starts at j = 2 because all the components of the first eigenvector (the one corresponding to the zero 
eigenvalue) are identical, which is the warranty of the final synchronization of the system. Here we can see 
the clear relation between topology (represented by the eigenvectors and eigenvalues of the Laplacian matrix) 
and dynamics. For long times all exponentials go to zero and the oscillators get synchronized. At short 
times the main contribution comes from small eigenvalues; then those nodes with similar projections on the 
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eigenvectors of small indices will get synchronized even at short times. In networks which are hierarchically 
organized this happens at all topological and dynamical scales. 

An additional significance of the importance of this relationship between spectral and dynamical proper- 
ties of identical oscillators comes from the work in fssj . The authors propose a method of network reduction 
based on the similarity of eigenvector projections. From the original network, nodes are merged according to 
the similarity of their components in the different eigenvectors, producing a reduced network; this merging 
procedure basically preserves the original eigenvalues of the Laplacian matrix in the new coarse grained 
Laplacian matrix. The authors determine the best clustering of the nodes and show that the evolution of 
identical Kuramoto oscillators in the original network and according to the original Laplacian is equivalent 
to the evolution of the reduced network in terms of the reduced Laplacian matrix. 

The above results refer to situations in which networks have clearly defined community structure. The 
approach we have shown enables one to deal with different time and topological scales. In the current 
literature about community detection fiolfl^. the main goal is to maximize the modularity, see Eq. In 
this case the different algorithms try to find the best partition of a network. Using a dynamical procedure, 
however, we are able to devise all partitions at different scales. In [s^ it is found that the partition 
with the largest modularity turns out to be the one for which the system is more stable, if the networks 
are homogeneous in degree^ If the networks have hubs, these more connected nodes need more time to 
synchronize with their neighbors and tend to form communities by themselves. This is in contradiction with 
the optimization of the modularity that punishes single node communities. From this result we can conclude 
that the modularity is a good measure for community partitioning. But when dealing with dynamical 
evolution in complex networks other related functions different from modularity are needed. 

For real networks, it has been shown that the same phenomenology applies 54]. These authors studied a 
system of Kuramoto oscillators, Eq. (jlOp with aij — cr/ki, arranged on the nodes of two real networks with 
community structures, the yeast protein interaction network and the Autonomous System representation of 
the Internet map. Both networks have a modular structure, but differ in the way communities are assembled 
together. In the former one, the modules are connected diversely (as for the synthetic networks analyzed 
before), while in the latter one different communities are interwoven mainly through a single module. The 
authors found that the transition to synchrony depends on the type of intermodular connections such that 
communities can mutually or independently synchronize. 

Modular networks are found in nature and they are commonly the result of a growth process. Nev- 
ertheless, these structural properties can also emerge as an adaptive mechanism generated by dynamical 
processes taking place in the existing network, and synchronization could be one of them. In particular, in 
[60| the authors studied the evolution of a network of Kuramoto oscillators. For a coupling strength below its 
critical value, the network is rewired by replacing links between neighbors with a large frequency difference 
with links between units with a small frequency difference. In this case, the network dynamically evolves 
to configurations that increase the order parameter. Along this evolution they noticed the appearance of 
synchronized groups (communities) that make the structure of the network to be more complex than the 
random starting one. 

Very recently [6l|, a slightly different model was considered, where the dynamics of each node is governed 

by 

being bij the betweenness centrality of the link, F^ the set of nodes that are connected to i, and a a time- 
dependent exponent. The authors use this dynamical evolution to identify communities. The betweenness is 
used as a measure of community coparticipation, since links between nodes that are in the same community 



have low betweenness 2l|. Starting from a synchronized state, a is decreased from zero and then the 



corresponding interaction strength on those links is increasingly enhanced. An additional mechanism that 
adjusts frequencies between neighboring nodes causes the final state to show partial synchronization among 
nodes that are in the same community. 



^Here stability (relative) of a given structure is understood to be the ratio between the final and initial times a partition 
remains synchronized. In terms of the number of connected components in Fig. |4]it corresponds to the length of the plateaus. 
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3.1.6. Synchronization by pacemakers 

It is worth mentioning the existence of other different approaches to the synchronization of populations of 
Kuranioto oscillators. So far we have referred to populations where the oscillators are nearly identical in the 
sense that they can have slightly different frequencies. Whenever there is a subset of units that play a special 
role, in the sense that they have substantially different frequencies than the rest in the population or they 
affect some units but are not affected by any of them, one usually refers to them as pacemakers. The effect of 
pacemakers has been studied in regular networks, as for instance in one-dimensional rings, two-dimensional 
tori and Cayley trees [6^ . So far, the only approach in a complex topology has been performed in (63j . 
There, the authors considered a system of identical units (same frequency) and a singular pacemaker. For 
an ER network they found that for a large coupling the pacemaker entrains the whole system (all units with 
the same effective frequency, that of the pacemaker), but the phase distribution is hierarchically organized. 
Units at the same downward distance from the pacemaker form shells of common phases. As the coupling 
strength is decreased the entrainment breaks down at a value that depends exponentially on the depth of 
the network. This result also holds for complex networks, as for instance in WS or SF networks, although 
the analytical explanation is only valid for ER networks. 

3.2. Pulse-coupled models 

In parallel to the studies described so far, some other approaches to synchronization in networks have 
invoked models where the interaction between units takes the form of a pulse. In particular, much attention 
has been devoted to models akin to reproduce the dynamics of neurons, e.g. integrate-and-fire oscillators 
(IFOs). The basics of an IFO system is as follows. The phase dynamics of any oscillator i is linear in time 
d(f>i{t)/dt = 1 in absence of external perturbations. However, when the oscillator i reaches the threshold 
(l)i{t) = 1 it sends a signal (or pulse) to the rest of the oscillators to which it is connected, and relaxes to 
4>i{t) — 0. The pulse can be considered to propagate instantaneously or with a certain time delay r, and 
when it reaches other oscillators induces a phase jump (pj -\- A(0j). The effects of the topology on the 

synchronization phenomena emerging in a network of IFOs are at least as rich as those presented in Sect. 
13. H although far more difficult to be revealed analytically. The main problem here is that the dynamics 
presents discontinuities in the variable states that are difficult to deal with. Nevertheless, many insights 
are recovered from direct simulations and clever mappings of the system. From direct simulations the 
first insights pointed directly to certain scaling relations between the synchronization time and topological 
parameters of networks. In ER networks, the scaling relation between the time to needed to achieve complete 
synchronization T, the number of nodes A'^, and the number of links M, was found to be 



with a — 1.30(5) and (3 — 1.50(5). Comparing this synchronization process with the same system on a 
regular square lattice, one realizes that the time needed to synchronize a random network is larger, specially 
in sparse networks [64| . In between of these two extremal topologies, some WS networks with a rewiring 
probability p were studied and were found to expand the synchronization time more than the original regular 
lattice. However, it was first pointed out that an appropriate normalization of the pulses received by each 
neuron, rescales the time to very short values. This phenomenon of normalization of the total input signal 
received by each oscillator has been repeatedly used to homogenize the dynamics in heterogeneous substrates. 

IFOs in SW networks were revisited later in [6^ to study the possibility of self-sustained activity induced 
by the topology itself. Considering a unidirectional ring of IFOs with density p of random long-range directed 
connections, the authors showed that periodical patterns persist at low values of p, while long-transients of 
disordered activity patterns are observed for high values of p. Responsible for this behavior is a tradeoff 
between the average path length and the speed of activity propagation. For low p configurations, the 
distances in the networks decrease logarithmically with size, while the superposition of activities is almost 
the same than in the regular configuration, i.e. the same activity occurs but in a "smaller" network able to 
self-sustain its excitation. However for large p, the superposition of activity between excited domains plays 




(42) 
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also an important role, and then both effects make the synchronized self-sustained activity collapse, leading 
to disordered patterns. 



In 66] , networks of nonidentical Hodgkin-Huxley [67J elements coupled by excitatory synapses in random, 
regular, and SW topologies, were investigated for the first time. The parameters of the model neurons were 
kept to stay below the bifurcation point, until the input arrives and forces the system to undergo a saddle- 
node bifurcation on a limit cycle. The dynamics of the system ends up into a coherent oscillation or in 
the activation of asynchronous states. In absence of a detailed analysis of the mechanism that generates 
coherence, the simulations showed several effects of the topology on the dynamics, the most interesting of 
which is that achieving synchronization in regular networks takes longer compared to SW, where the existence 
of short-cuts favors faster synchronization. The results obtained in all cases show that the randomness of 
the topology has strong effects on the dynamics of these models, in particular the average connectivity is a 
control parameter for the transition between asynchronous and synchronous states. In Fig. [5] we present a 
phase diagram for the Hodgkin-Huxley model in SW networks with varying average degree (k) and rewiring 
probability p. A detailed analysis of sparse random networks of general IFOs was exposed in [H^ . Their 
analytic results are in agreement with the previous observations. Very recently [69|, a SW network of 
non-identical Hodgkin-Huxley units in which some of the couplings could be negative was analyzed; they 
surprisingly found that a small fraction of such phase-repulsive links can enhance synchronization. 

In a slightly different scenario 70] a system of pulse-coupled Bonhoeffer-van der Pol-FitzHugh-Nagumo 
oscillators [71] in WS networks was studied numerically. This study reports a major influence of the average 
path length of the network on the degree of synchronization, whereas local properties characterized by 
clustering and loop coefficients seem to play a minor role. In any case, the authors warn that the results are 
far from being conclusive, since single characteristics of the network are not easily isolated. We will come 
back to this issue in the next section, when dealing with the stability of the synchronized state. 

The works reviewed so far in this subsection are based on the assumption that the coupling is fixed, 
and that the only source of topological complexity is embedded in the connectivity matrix. The authors in 
[73 ] showed that for networks of pulse-coupled oscillators with complex connectivity, coupling heterogeneity 
induces periodic firing patterns, which replace the state of global synchrony. The coupling heterogeneity has 
a critical value from which the periodic firing patterns become asynchronous aperiodic states. These results 
are in agreement with the observations described in previous works and allow us to state that a certain 
degree of complexity in the interaction between pulse coupled oscillators is needed to observe regular (or 
ordered) patterns. However, once a critical level of complexity is surpassed, asynchronous aperiodic states 
dominate the dynamic phenomena. 

3.3. Coupled maps 

Maps represent simple realizations of dynamical systems exhibiting chaotic behavior. At a flrst sight 
they can represent discrete versions of continuous oscillators. Coupled populations of such rather simple 
dynamical systems have been one of the paradigmatic models to explain the emergence and self-organization 
in complex systems due to the rich variety of global qualitative behavior they give rise to. From a more 
practical point of view coupled map systems have found a widespread range of applications, ranging from 
fluid dynamics and turbulence to stock markets or ecological systems [l| . Since these systems are nowadays 
known to have complex topologies, populations of maps coupled through a complex pattern of interactions 
are natural candidates to study the onset of synchronization as an overall characteristic of the population. 

Coupled maps have been widely analyzed in regular lattices, trees and also in glob al connectivity schemes. 
The first attempt to consider connectivities in between these extreme cases is |73|. He proposed a system 
formed by units, whose individual dynamics are given by the logistic map, that are connected to a flxed 
number k of other units randomly chosen (multiple and self-links are permitted). The evolution rule for the 
units is 



i 

A linear stability analysis of this system is performed in terms of the eigenvalues of the matrix A. For the 
logistic map it is shown that for fc > 4 the maps synchronize. The time the system needs to synchronize 
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Figure 5: Phase diagram which shows the regions of oscillatory (clear) and nonoscillatory (dark) activity of the network in the 
(fc,p) plane. The island that appears on the right side indicates that the SW (for some range of values of k) is the only regime 
capable to produce fast coherent oscillations in the average activity after the presentation of the stimulus. Prom [66( . 

decreases with the connectivity k and also with the system size, ahhough in the latter case the time saturates 
for large values of the system size. When the connectivity pattern is changed to a modified WS model (by 
adding long-range short-cuts but not rewiring), the authors in [tJ] showed that just a nonzero value of the 
addition probability is enough to guarantee synchronization in the thermodynamic limit. 

In another early attempt to include non-regular topologies in chaotic dynamics [tsI , a SW network was 
analyzed, in which 

N 

e.,it + !) = (!- a)f {em + E«»^/ (^^W) ' (44) 

where k is the number of shortcuts in the network, and a the coupling constant. Each unit evolves according 
to a sine-circle map |76j 

f{e) = e + —Bin{2'Ke) (modi), (45) 

which provides a simple example for describing the dynamics of a phase oscillator perturbed by a time- 
periodic force. Here X is a constant related to the external force amplitude and < O < 1 is the ratio 
between the natural oscillator frequency and the forcing frequency. It is observed that synchronization, 
in terms of a parameter related to the winding number dispersion, is induced by long-range coupling in a 
system that, in the absence of the shortcuts, does not synchronize. 

A slightly different approach was conducted in [77[, who considered a population of units evolving 
according to 

x,{t -I- 1) = (1 - a)f {x.,{t)) + YT.f (^^- W) ■ (46) 

' jer. 

They obtain the stability condition of the synchronized state in terms of the eigenvalues of the normalized 
Laplacian matrix [Sij — aij/ki) and the Lyapunov exponent of the map f{x). Furthermore, they also find a 
sufficient condition for the system to synchronize independently of the initial conditions, namely 

(l-aA2)sup|/'| < 1, (47) 

where A2 is the smallest non-zero eigenvalue of the normalized Laplacian matrix. They demonstrate their 
results for regular connectivity patterns as global coupling and one-dimensional rings with a varying number 
of nearest neighbors, since in these cases the eigenvalues can be computed analytically. Complexity in the 
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connectivity pattern is introduced in different ways. In these cases one needs numerical estimates of the 
eigenvalues to compare the synchronization condition (j47p with the simulation of the model (|46p . By using 
a quadratic map f{x) = 1 — ax'^ [t^ and choosing the free parameter a in a range where different regimes 
are realized, they find that in a random network the system synchronizes for an arbitrary large number of 
units, whenever the number of neighbors is larger than some threshold determined by the maximal Lyapunov 
exponent. This implies a remarkable difference to the one-dimensional case where synchronization is not 
possible when the number of units is large enough. For a WS model their main finding is that a quite high 
value of the rewiring probability [p > 0.8) is needed to achieve complete synchronization. Finally for BA 
networks the behavior is comparable to the random ER case. 

Following a similar line, in [t^ it is studied the behavior of a model where the interaction between the 
units can be strengthened according to the degree. In this case 

xt+i,, = (1 - a)f{xt,^) + E ''ffM, (48) 

where A/i — X^jer appropriate normalizing constant. Here the function f{x) is also a quadratic 

map. The authors study first BA networks. When a = the model is equivalent to that discussed previously; 
in this case they find the existence of a first-order transition between the coherent and the noncoherent 
phases that depends on both the mean connectivity and the coupling a. As varying a, the parameter of the 
quadratic map, they find that these two critical values are related by the power law CTc oc k~'^. The effect of 
a being larger than zero is only quantitative, since in this case the transition appears at smaller values of the 
interaction as compared with the usual case. Additionally, the authors studied two types of deterministic SF 
networks: a pseudofractal SF network introduced in and the Apollonian network introduced in |80(. In 
both cases, there is no coherence when a = and a — 2. This fact leads the authors to conclude that some 
degree of randomness that shortens the mean distance between units is needed for achieving a synchronized 
state, since in these networks the SF nature is not related to a short mean distance. Nevertheless, this 
situation is avoided if the contributions from the hubs are strengthened, by making a > 0. 



Another set of papers deals with units that are coupled with some transmission delay [81|, l82, l83( . For 
instance, in ^81.] the authors propose a model in which all units have the same time delay (in discrete units) 
with respect to the unit considered: 

x.it + 1) = / (.x,(t)) + [/(^^- - ^^i)) - f (^^w)] ■ (49) 

For a uniform delay Ty = r Vz,j, they show analytically, and numerically, that the delay facilitates syn- 



chronization for general topologies. In any case, this fact confirms the results obtained in [77| that ER and 
SF networks are easier to synchronize than regular or SW ones. Furthermore, one of the implications of 
connection delays is the possibility of the emergence of new collective phenomena. In 82, 8^ the authors 



considered uniform distributions of (discrete) time delays. Their main result is that in the presence of ran- 
dom delays the synchronization depends mainly on the average number of links per node, whereas for fixed 
delays there is also a dependence on other topological characteristics. 



In a more general framework [8J, |85|, |86| , the following problem is considered 

x,{t + !) = (!- a)f {x,{t)) + E 5 i^jit)) ■ (50) 

' jer. 

For a logistic map g{x) — fix{l — x) (although analyses on other maps as the circle map and the 
tent map have also been performed) the authors show a phase diagram in which the different stationary 
configurations are obtained as a function of the coupling strength a and the parameter of the map fi. 
The stationary configurations are classified in the following way: turbulence (all units behave chaotically), 
partially ordered states (few synchronized clusters with some isolated nodes), ordered states (two or more 
synchronized clusters with no isolated nodes), coherent states (nodes form a single synchronized cluster), 
and variable states (nodes form different states depending strongly on initial conditions). The critical value 
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Figure 6: (color online) Left: Distribution of the topological distances for a tree with 1000 nodes (bullets) and the distance 
inside the synchronized clusters (other symbols) for a = 0.017 and different initial conditions. Right: visualization of the tree 
with five interconnected clusters of synchronized nodes marked by different colors. From [87| . 



of the coupling above where phase synchronized clusters are observed depends on the type of network and 
the coupling function. As a remarkable point, it is found that two different mechanisms of cluster formation 
(partial synchronization) can be distinguished: self-organized and driven clusters. In the first case, the nodes 
of a cluster get synchronized because of intracluster coupling. In the latter case, however, synchronization 
is due to intercluster coupling; now the nodes of one cluster are driven by those of the others. For a linear 
coupling function g(x) = x, self-organization of clusters dominates at weak coupling; when increasing the 
coupling strength, a transition to driven-type clusters, almost independent of the type of network, appears. 
However, for a nonlinear coupling function the driven type dominates for weak coupling, and only networks 
with a tree-like structure show some cluster formation for strong coupling. 

Finally, it is worth mentioning the very Ref. [s^l, in which the authors consider a SF tree (preferential- 
attachment growing network with one link per node) of two-dimensional standard maps: 

x' = X + y + /i sin(27rx) (modi) , ^ 

y' = y + fj,sin{2'Kx). ^ ' 

The nodes are coupled through the angle coordinate (x) so that the complete time-step of the node i is 

x^{t + l) = il-e)x'S) + i-EjerA^At)-^'M 

y^it + ^) = i^-e)ym- ^ ' 

Here, (') denotes the next iteration of the (uncoupled) standard map ([5T|) and t denotes the global discrete 
time. The update of each node is the sum of a contribution given by the updates of the nodes, the ' part, 
plus a coupling contribution given by the sum of differences, taking into account a delay in the coupling from 
the neighbors. By keeping fi = 0.9 such that the individual dynamics is in the strongly chaotic regime, the 
authors analyze the dependence on the interaction strength a. For small values of the coupling the motion 
of the individual units is still chaotic, but the trajectories are contained in a bounded region. With further 
increments of the coupling, the units follow periodic motions which are highly synchronized. In this case, 
however, synchronization takes place in clusters, each cluster having a common value of the band center 
around which the periodic motion occurs, and center values appear in a discrete set of possible values. These 
clusters form patterns of dynamical regularity affecting mainly nodes at distances 2, 3, and 4, as shown in 
Fig. [5{left). In fact, the histograms of distances between nodes along the tree and between nodes belonging 
to the same synchronized cluster have different statistical weights only for these values of the distances. 

All previously discussed chaotic models are discrete time maps, which are appropriate discrete versions 
of chaotic oscillators. Nevertheless, we also notice a couple of works dealing with time continuous maps. 
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In [88[ the authors analyze a WS network of Rossler osciUators, where the parameters are chosen to ensure 
that the system generates chaotic dynamics. The basic observation is that the network synchronizes when 
the couphng strength is increased, as expected. Another interesting result is that the mean phase difference 
among the chaotic oscillators decreases with the increasing of the probability of adding long-range random 
short-cuts. Along the same line, in [s^ it is considered a system of Rossler oscillators on BA networks. The 
only tuning parameter of the BA networks is to, the number of links that a newly added node has. For to = 1 
(SF trees) there is no synchronized state for a large number of oscillators. Increasing m synchronization is 
favored. The topological effect of increasing to is to create loops, but it is shown that this is not the only 
fact that improves synchronization. 

Finally, a different and interesting proposal was made in [o^l where a fixed 1-d connectivity pattern 
is complemented by a set of switching long-range connections. In this case, it is proven that interactions 
between nodes that are only sporadic and of short duration are very efficient for achieving synchronization. 

As a summary, we can say that most of the works deal with particular models of coupled maps (logistic 
maps, sine-circle maps, quadratic maps, ...). Thus, it is possible, to obtain in some cases not only the 
conditions of local stability of the completely synchronized state but the conditions for the synchronization 
independently of the initial conditions. In general, the addition of short-cuts to regular lattices improves 
synchronization. There are even some cases, for which synchronization is only attainable when a small 
fraction of randomness is add to the system. On the contrary, in the next section we will discuss the linear 
stability of the synchronized state for general dynamical systems. 



4. Stability of the synchronized state in complex networks 

In the previous section we have reviewed the synchronization of various types of oscillators on complex 
networks. Another line of research on synchronization in complex networks, developed in parallel to the 
studies of synchronization in networks of phase oscillators, is the investigation of th e stability of the com- 
jletely synchronized state of populations of identical oscillators. The seminal work by Barahona and Pecoral 



9ll | initiated this research line by analyzing the stability of synchronization in SW networks using the Master 
Stability Function (MSF). The framework of MSF was developed earlier for the study of synchronization 



of identical oscillators on regular or other simple network configurations [92, l93|. The extension of the 
framework to complex topologies is natural and important, because it relates the stability of the fully syn- 
chronized state to the spectral properties of the underlying structure. It provides with an objective criterion 
to characterize the stability of the global synchronization state, from now on called synchronizability of net- 
works independently of the particularities of the oscillators. Relevant insights about the structure-dynamics 
relationship has been obtained using this technique. 

In this section, we review the MSF formalism and the main results obtained so far. Note that the MSF 
approach assesses the linear stability of the completely synchronized state, which is a necessary, but not a 
sufficient condition for synchronization. 



4-1. Master Stability Function formalism 

To introduce the MSF formalism, we start with an arbitrary connected network of coupled oscillators. 
The assumption here for the stability analysis of synchronization is that all the oscillators are identical, 
represented by the state vector x in an m-dimensional space. The equation of motion is described by the 
general form 

x = F(x). (53) 

For simplicity, we consider time-continuous systems. However, the formalism applies also to time-discrete 
maps. We will also assume an identical output function H(x) for all the oscillators, which generates the 
signal from the state x and sends it to other oscillators in the networks. In this representation, H is a 
vector function of dimension to. For example, for the 3-dimensional system x = {x,y,z), we can take 
H(x) = (x, 0,0), which means that the oscillators are coupled only through the component x. H(x) can 
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be any linear or nonlinear mapping of the state vector x. The N oscillators, i = 1, . . . , iV, are coupled in a 
network specified by the adjacency matrix A ~ (ay )- We have 

N 

X, = r(xi) + cr^ayU;y[H(xj) - H(xi)] (54) 

N 

= F(x,) -a^GyH(xj), (55) 

being Wij > the connection weights, i.e., the network is, in general, weighted. The coupling matrix G is 
Gij = —aijWij lii ^ i and Gu — X^jLi O'ijWij- When the coupling strength is uniform for all the connections 
{wij = 1), the network is unweighted, and the coupling matrix G is just the usual Laplacian matrix L. By 
definition, the coupling matrix G has zero row-sum. Thus there exists a completely synchronized state in 
this network of identical oscillators, i.e., 

Xi(t)=X2(i) = ... = XAr(i) = s(i), (56) 

which is a solution of Eq. (j55p . In this synchronized state, s{t) also approaches the solution of Eq. (j53|l . i.e., 
s — F(s). This subspace in the state space of Eq. ([55]) . where all the oscillators evolve synchronously on the 
same solution of the isolated oscillator F, is called the synchronization manifold. 

4.1.1. Linear Stability and Master Stability Function 

When all the oscillators are initially set at the synchronization manifold, they will remain synchro- 
nized. Now the crucial question is whether the synchronization manifold is stable in the presence of small 
perturbations 5xi. To assess the stability, we need to know whether the perturbations grow or decay in 
time. The linear evolution of small Sxi can be obtained by setting Xi{t) = s{t) + Sxi{t) in Eq. ([55|) . 
and expanding the functions F and H to first order in a Taylor series, i.e., F(xi) — F(s) + DF{s)Sxi and 
H(xi) = H(s) + Z)H(s)(5xi. Here DF{s) and Z3H(s) are the Jacobian matrices of F and H on s, respectively. 
This expansion results in the following linear variational equations for Sxi, 

N 

6±^ = DF(s)(5x, - (tDU{s) ^ Gi^Sxj. (57) 

j=i 

The variational equations display a block form, each block (ij) having m components. The main idea here is 
to project Sx into the eigenspace spanned by the eigenvectors of the coupling matrix G. This projection 
can operate in block form without affecting the structure inside the blocks. By doing so, Eqs. ([57]) can be 
diagonalized into N decoupled eigenmodes in the block form 

ii^[DFis)-aXiDll{s)]^i, l = l,---,N, (58) 

where is the eigenmode associated with the eigenvalue A; of G. Since G has the property of zero row-sum, 
the minimal eigenvalue is always zero, i.e., Ai = 0, with the corresponding eigenvector Vi = (1, 1, . . . , 1). So 
the first eigenmode ^1 — £'F(s)^i corresponds to the perturbation parallel to the synchronization manifold. 
The other — 1 eigenmodes are transverse to the synchronization manifold and should be damped out to 
have a stable synchronization manifold. 

In general, the eigenvalues spectrum is complex when the coupling matrix is not symmetrical, e.g., 
when the network is directed (a.y 7^ aji) or when the coupling weights are asymmetrical (wij 7^ Wji). 
In the literature, the characterization of network topology and the analysis of synchronization and other 
dynamical processes have mainly focused on undirected and unweighted networks. This case corresponds to 
a symmetric G and therefore all its eigenvalues are real, which makes the analysis simpler. In the following 
we also assume that all the eigenvalues are real, which is always the case for symmetric G but it can also 
be true for non-symmetric cases, as will be discussed later. 
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When the eigenvalues spectrum is real, it has the following properties: (i) Ai = due to zero row sum 
of G; (ii) A; > since G is positive semidefinite and (iii) there is only one zero eigenvalue if the network is 
connected. Accordingly, the eigenvalues can be ordered as 

= Ai < Aa • •• < Aat. (59) 

In a connected network of N oscillators, A2 and Aat are then, the minimal and the maximal non-zero 
eigenvalues, respectively. 

Importantly, we observe that all the individual variational equations in the system of equations Eq. (I58|) 
have the same form 

e = [DF{s) - aL»H(s)] ^ (60) 

They only differ by the parameter ai = aXi. Now if we know the stability of the solution ^ = for any 
reasonable value of a, then we can infer the stability for any eigenmode with ai = aXi. To assess the stability 
of this master variational equation (|60p . we calculate its largest Lyapunov exponent Amax as a function of a, 
the resulting function is the master stability function. The evolution of small ^ is then described on average 
as ~ exp[Amax(Q;)i], and the mode is stable with ||^|| if A„iax(a) < 0. 

So far, we have assumed that the eigenvalues A; are real, and the problem is to compute the MSF Aniax(a) 
for real values of a. For the general case of complex eigenvalues A; one needs to evaluate the MSF in the 
complex plane a = {aji,ai). This scenario is mathematically more intricate and less results are available. 



For many oscillators types Amax(a) < in a bounded region of this plane j92l. l93l|. However, semi-bounded 
(generally speaking unbounded) regions where Amax(a) < can also occur for some specific oscillators. We 
will see later how this difference between region bounds affects the measurement of synchronizability. 

The reader may find the framework of master stability function abstract. Here we present it in a 
physically intuitive manner with the help of a schematic diagram. For this, we restrict ourself to the case 
where a is real. Let us consider first two coupled oscillators. The first one is autonomous and evolving along 
the trajectory s(t), and the second one is driven by the first one with a coupling strength a, as depicted in 
Fig. [Tl^a). The dynamical equations read, 

s = F(s), (61) 
x = F(x) +a[H(s) - H(x)]. (62) 

Then immediately, we obtain the linear variational equation for the synchronization difference ^(t) = x(i) — 
s{t) as in Eq. (|5D|) . This means that the MSF describes the stability of the state in which the two oscillators 
are synchronized. 

In a similar spirit, the equations for the — 1 transverse eigenmodes in Eq. (|58p can be graphically 
represented as iV — 1 oscillators driven by a common autonomous oscillator s = F(s), and the corresponding 
coupling strengths are aXi, as depicted in Fig. [71(b). In this representation, the mode decomposition decouples 
the complex network connections into many pairs of oscillators and the stability thus can be understood 
from that of the two coupled oscillators, the MSF. The complete synchronization state s is stable when all 
the N — 1 oscillators in Fig. [Tljb) are synchronized by the common forcing signal s. This picture will be 
useful later on in this section. 



4- 1.2. Measures of synchronizability 

In the case of real eigenvalue spectrum of G, the MSF only presents real values of a. Here we distinguish 
two cases that will lead to different criteria for synchronization in networks. The first case refers to bounded 
MSFs, where Ainax(a) < within a finite interval ai < a < a2 (Fig. [S]). A physical picture of this case for 
two coupled oscillators, as those in Fig.[7Ja), is that the driven oscillator will be synchronized, x(i) — s{t), 
if the coupling strength a S (01,0:2). Bounded MSFs are quite common for many oscillators F and many 
coupling functions H. 

For discrete time chaotic maps, the MSFs are always bounded, the reason is the following: in discrete 
time systems, Eq. ((50)) reads as 



= [-DF(s(n)) -ai?H(s(n))]e„ = J(a,s(n))C„ = J"^o, 
27 
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Figure 7: Schematic plot of the eigenmode decomposition, (a) For 2 unidircctionally coupled nodes; (b) for a network of TV 
coupled nodes. 




Figure 8: Four master stability functions for coupled Rossler oscillators: chaotic (bold) and periodic (regular lines); with y 
coupling (dashed) and x coupling (solid lines). (The curves are scaled for clearer visualization). From |91| . 
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where n is the iteration step and J" = n"^g J(a, s(/)). The largest Lyapunov exponent is [94| 



Amax= lim iln^= hm ;^ln||J^||. 

T^ooT 1 1^0 1 1 T^oaT 

where T is the number of iterations. Now suppose a > is a very large value, then at each step n, 
J{a,s{n)) sa —aD'H.{s{n)), and J'^ ~ (— a)-^n^^Q £'H(s(n)). As a result, A^ax ~ Ina + ^. Here 
fj, = limr^oo In 1 1 1 1 is a finite measure on the chaotic attractor similar to the largest Lyapunov ex- 
ponent Af of the isolated chaotic map F, Af — limT^oo y In || Jj||, where — n^~Q£>H(s(n)) and 
Jp — n^ZQDF(s(n)), respectively. Consequently, Amax > when a is large enough, i.e. the MSF always 
becomes positive when the coupling strength is large enough, thus it is bounded. 

There are also situations where the MSFs are unbounded, and Amax < for a > ai without the upper 
limit a2- This happens, for example, in time-continuous systems with H(x) — x where the oscillators are 
linearly coupled through all the m corresponding components. In this case, the driven oscillator in Fig.[71[a) 
will be synchronized if the coupling strength a > ai. 

To synchronize a network of oscillators, all the iV — 1 oscillators in Fig. [T^b) must be synchronized by the 
common forcing signal s{t). In the case of bounded MSFs, this requires that aXi € (ai, (22) for 2 < ^ < A^. 
Explicitly, the following condition is necessary for the stability of the synchronization state of the network, 

ai < cr A2 < (tAs < • • • < (tAat < a2. (63) 

This condition can be only fulfilled, for some values of a, when the eigenratio R satisfies the following 
relation 

R^^<^. (64) 
A2 ai 

Therefore, we conclude that it is impossible to synchronize the network if i? > a2/ai, since there is no a for 
which the fully synchronized state is linearly stable. On the contrary, if i? < a2/ai the synchronous state 
is stable for CTmin < <J < (Tmax where (Jnim = ai/A2 and cr,nax = 0^2 /Aat, respectively. 

If the MSF is unbounded, then the synchronization manifold is stable if it exists a satisfying 

ai < 0X2 < crAa < • • • < (jXn, (65) 

which will be true for any a larger than a certain synchronization threshold (Xmin- Thus, the larger A2 the 
smaller synchronization threshold (Tmin Moreover, as previously discussed in Sect. 13.1.41 A2 also 



plays a special role in the time needed to achieve complete synchronization |52l |. 

Note that the eigenratio R and the nonzero minimal eigenvalue A2 depend only on the network structure, 
as defined by G. For bounded MSF, if R is small, the condition in Eq. (|64p is, in general, easier to satisfy. 
From this, it follows that the smaller the eigenratio R the more synchronizable the network and vice versa 



9l( . In this sense, we can characterize the synchronizability of the networks with R and A2, without referring 
to specific oscillators. As will be discussed soon, these two types of synchronizability on a specific network 
can depend on different parameters and then can be correlated to different network descriptors. In this 
review, synchronization based on the eigenratio R is called Type I and synchronizability based solely on A2 
is called Type 11 synchronizability. 

For general directed networks, the spectrum of the coupling matrix G is complex, and it is still unclear 
how to develop simple measures, without referring to the MSF of specific oscillators and coupling functions, 
to evaluate the synchronizability of different (directed) networks. This is probably one of the reasons why 
the vast majority of works using the MSF have focused on undirected and unweighted networks whose 



spectra are real. In a recent work [97[, both the ratios of the real and imaginary parts of the eigenvalues 
are used. However, it should be noted that networks with the same ratios (or even the same minimal and 
maximal real and complex parts), but different boundaries of eigenvalues in the complex plane, will have 
different synchronization thresholds for the same MSF. This is contrast to the cases of real spectra, where 
the synchronization thresholds for the same MSF will be identical for networks with the same A2 and Aat, 
irrespective to the other spectral properties. 
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Figure 9: (color online) Decay of the eigenratio i? in a A'^ = 100 lattice as a fraction / of the N{N — l)/2 possible edges are 
added following purely deterministic, semirandom, and purely random schemes. Networks become synchronizable below the 
dashed line (/3 = 012 1 a\, where 01 and 02 are from Fig.[8ll. The squares (numerical) and the solid line [analytic Eq. II66I 1] show 
the eigenratio decay of regular networks through the deterministic addition of short-range connections. The dot-dashed line 
corresponds to purely random graphs [Eq. ([TTJ], which become almost surely disconnected and unsynchronizable at / ~ 0.0843. 
The semirandom approach (dots, shown for ranges fc = 1,2,4,6,10,14) is more efficient in producing synchronization when 
k < InN. From 



The key advantage of the MSF framework is that it provides an objective criteria ( A2 and Aat) to assess 
the synchronizabihty of complex networks without referring to specific oscillators. The drawback is that it 
only informs about the dynamics towards synchronization from small perturbations of the synchronization 
manifold (linear stability). The study of synchronizabihty is then converted into the investigation of the 
eigenvalues of the coupling matrix G of the networks. In the current scenario of the MSF, natural questions 
are: What is the synchronizabihty of different types of complex networks? Which structural properties are 
related to or control the synchronizabihty? 

Note that although the more appropriate approach seems to be that of following the general framework 
of graph theory to investigate the spectral properties of networks, many authors have tried to relate a 
single statistical property of networks with synchronizabihty, sometimes misinterpreting and generalizing 
the results without a proper estimation of their constraints. In the following section, wc will summarize 
results on the synchronizabihty of typical network models and then describe the relationship between graph 
theoretical measures and the synchronizabihty of complex networks. Further insights have been obtained 
in the scope of graph theory by providing with bounds that constrain the eigenvalues of the Laplacian of 
networks. 



4-.1.3. Synchronizability of typical network models 

We will present the main findings related to synchronizability in the most common classes of complex 
networks found in the literature: Regular, SW, ER, and SF networks. 



Regular networks. For a long time, synchronization of chaotic oscillators has been studied on regular net- 
works 92, 9^. A typical example of a regular network is a cycle (or ring) of N nodes each coup led to its 2k 
nearest neighbors with a total of Nk links. The eigenvalues of the coupling matrix are 9^, 9l| 



A, = 2fc 



A , 2^(/-l)j\ 

2z.^°n — ^ — , 



(66) 
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By a series expansion, we obtain 



A.. ^''^"-;;P+" .A„ = (2t+l)(l+2/3.). (67) 



and the eigenratio R for fc <C -/V can be approximated as 

These relations show that regular networks have a poor Type I and Type II synchronizability. The asymp- 
totics on the system size A2 ^ 1/^^, makes synchronization effectively impossible for both bounded and 
unbounded MSFs of chaotic oscillators with ai > 0. In such regular networks, usually we have complicated 
pattern formation (waves) instead of a stable spatially homogeneous (completely synchronized) state [99|. 
A way to improve this situation is, for fixed iV, adding connections to generate a higher range fc, since A2 
increases approximately as A2 ~ the the eigenratio R decreases as i? ~ for k <^ N, and synchroniz- 
ability of Type I and II is enhanced. 

SW networks. In joij , SW networks are obtained by adding NS links at random, to a regular network where 
each node is connected to its 2k nearest neighbors, so that the average number of shortcuts per node is S. 
As shown in Fig. [9l adding a small fraction of such random connections reduces the eigenratio R so that the 
synchronizability is improved significantly. 

The eigenvalues of the SW networks have been obtained through a perturbation analysis of the SW 
Laplacian L — ^U' . Here L° is the deterministic Laplacian of the regular networks and U' the Laplacian 
formed by the random shortcuts. In the stochastic Laplacian matrix (symmetric, zero row-sum), any of 
the remaining entries N{N — 2fc — 1) of L'^ takes the value 1 with probability ps = 2S/ {N — 2fc — 1) and the 
value with probability (1 — Ps)- For ps ^ 1, and N-^^^ < k N, the perturbations of the eigenvalues are 

£A^'^ ~ Nps - v/37rp,/4, eA^^ ~ Np^ + v/37rp,/4. (69) 
The extreme eigenvalues are 

A2 = a(")+£AW, Ajv = Ai^)+eAi^\ (70) 

where Aj^'* and A^-* are the eigenvalues of the regular networks (L") as in Eq. (p7|) . From this analysis, it 
follows that for a fixed small value of S, the minimal non-zero eigenvalue A2 is driven away from X^2^ « 
to eX'i^ « Np, « 2S for any SW network with large N and N-^^^ < k <$:. N, while the maximal eigenvalue 
Xn is not affected very much for small S. This means that both types of synchronizability are mainly 
determined by the average number of shortcuts per node S. 

In SW networks, the va rianc e of the degree distribution raises as the regular network is rewired with 
an increasing probability p jlOd ] or when more shortcuts are added, see Fig. fTUT b). This process results 
in an improvement of the synchronizability (the eigenvalue ratio, R, is reduced), as illustrated in Fig. [9] 
and Fig. fTOF a). This is because A2 increases proportionally to the number of shortcuts per node, i.e., 
A2 ~ 25 = 2kp. 

Random networks. In purely random graphs, in which a fraction / of the N{N — l)/2 possible links is 
established at random, the eigenratio is a function of / and N [9l[. It reads 



^ ^ Nf + ^2/(1 - fjm^ .7^. 

Nf^^2fil-f)N\nN' 

Note that the networks are synchronizable only when / > 2\nN/{N + 21niV), where it is almost sure 
that the networks are connected [lOlj. If this condition is verified the synchronizability is improved with 
increasing /, as seen in Fig. [S) 
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Figure 10: Synchronizability (eigenratio R, (a)) and heterogeneity of degrees (variance (t| of the degree distribution (b)) of SW 
networks as a function of th e rew iring probability p. Inset of (b): average path length £ vs p. The network has a size = 2000 
and the range k = 3. From [lOO( |. 




Figure 11: Synchronizability of random SF networks. Plotted are eigenvalues A2 (circles) and Ajv (squares) and eigenratio R 
(triangles), (a) As functions of the exponent 7 for network size N = 2^"; open symbols for fc^in = 5 and filled symbols for 
kmin = 10. (b) As a functions of network size A'^ at 7 = 3 and fcmin = 10. The results are averaged over 50 realizations. 



Note that for small / < \nN/N, the purely random networks are almost surely disconnected and thus 
non-synchronizable. On the contrary, the regular backbone of nearest connections with fc > 1 can already 
make the semirandom SW networks connected regardless of N and thus synchronizable as a whole. In this 
sense, semirandom SW networks turn out to be much superior to the purely random networks in terms of 
synchronizability, as can be seen from Fig. [9] for small fc, fc < In (fc = 1 ^ 4 in networks with N = 100). 
The improvement is even more pronounced for larger N. 

SF networks. Here we discuss the synchronizability of SF networks for diffe rent values of the degree distri- 



bution exponent 7. To this end, we consider a random model introduced in 102l | to construct SF networks. 
The algorithm works as follows, first, a degree ki is assigned to each node i according to the probability 
distribution P(fc) ^ k~^ and fc > fcmin- The network is generated by randomly connecting the nodes so that 
node i has exactly the prescribed ki links to other nodes, prohibiting self- and repeated links. 

The dependence of the eigenvalues A2 and Xn and the eigenratio i? on 7 and N are shown in Fig. [TlJ 
We observe that A2 has no noticeable dependence on 7 and N. However, Xn becomes larger as the degree 
heterogeneity is increased. So the changes of the eigenratio R follow the trend of Ajv closely. This result 
is somehow expected given that the largest eigenvalue A at is intimately related to the degree of hubs, and 
this is the essential fingerprint of SF networks with different exponents 7. On the other hand, A2 increases 
with fcmin, and R is larger at smaller fcmin for the same 7 and N. The variation of i? as a function of 7 
was reported in |l03l |. The dependence of the synchronizability on fcmin and N in this random SF network 
model turns out to be very similar in the BA growing model, as shown in [l04]. The conclusion is that in SF 
networks, the two types of synchronizability (I and H) associated to R (bounded MSF) and A2 (unbounded 
MSF) are very differen t. 

In an early paper |l05l | that studied robustness and fragility of synchronization of SF networks, it was 
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reported that A2 is a constant almost unrelated with fcmin (^min = 3,5,7,9,11). This result in [105| is 
inc onsi stent with the observation in Fig. 1111 with the work in 104], and with graph theoretical analysis 
in |l06l |. which we will discuss in more detail later on. Thus the observation of robustness and fragility of 
synchronizability in [105] (changes of A2 due to random or deliberate attack of the nodes) should be taken 
cautiously, and a detailed reexamination of this issue is mandatory before assessing conclusions. 

4-. 1.4- Synchronizability and structural characteristics of networks 

The relationship between structural characteristics of networks and synchronizability has been explored 
intensively in the literature, mainly based on numerical experiments on various network models. The 
observations, which are summarized in what follows, are quite confusing. The main problem is that many 
works have made a naive use of complex network models to assess synchronizability. Very often network 
models do not allow us to isolate one structural characteristic while keeping the other properties fixed. For 
this reason many results have been misinterpreted. As we will see in the next section, a more objective graph 
theoretical analysis sheds some light to the whole problem. Unfortunately, the complete understanding of 
the structure-synchronizability relationship is still missing. 

Synchronizability dependence on £. The average shortest path length £ is a property of the network closely 
related to the efficiency of information processing. Most real-world complex networks are characterized by 
a small ^ < In [s^]. Indeed, it has been conj ectured and rationalized that in biological neuronal networks, 
£ has been minimized by evolution llOTL llOSi] . Generally speaking, i is lower in SF networks than in ER 
networks due to the presence of hubs [109j , and £ is lower in SW networks than in regular lattices due to 
the presence of shortcuts. 

In [§] it is suggested that the decrease in the distance in the WS network would lead to more efficient 
coupling and thus enhanced synchronization of the oscillators. Investigation of phase oscillators [s^ ] or 
circle maps [75] on WS networks has shown that when more and more shortcuts are created at larger 
rewiring probability p, the transition to the synchronization regime becomes easier. On the other hand, 
the synchronizability of identical oscillators follows the same trend of £, in networks with fixed N and k, 
as p is increased (Fig. [TU)) . A similar type of behavior is observed if shortcuts are added to the regular 
networks (Fig. H]). From these observations, the generalized conclusion that smaller distances will always be 
correlated to enhanced synchronization, has been intuitively used by many authors, and is not so. Indeed, 
a more detailed analyses of various network models have shown that there is no direct relationship between 
£ and the synchronizability of the networks. The reason is that the transition to the small- world regime 
occurs at a value of the rewiring probability for which there is no significant effect on A2. 

In fact, in WS networks, is a function of the ne twor k size N, the degree of the nodes in the original 
regular network, fc, and the randomness parameter p |ll0l ] 

£{N,k,p)^jf{pkN), (72) 

where f{u) is a universal scaling function, f{u) = const if u ^ 1 and f{u) — ln(u)/u if m ^ 1. From 
this result, it turns out that £ begins to decrease with p, and consequently the SW behavior emerges, for 
P ?i Psw = l/Nk. At p = psw the average number of shortcuts per node is 5 ~ 1/A^, and then A2 ~ 1/A^ as 
well. This shows that at this point the synchronizability is not enhanced by the rewiring. To achieve such 
an enhancement, the density of shortcuts has to be independent of N, which happens for p > Psync = l/fc, 
that is deep in the SW regime. In other words, in the intermediate region pgw < p < Psync, ^ decreases 
wh ile the synchronizability o f the system remains roughly the same. 



Barahona and Pecoral [91| showed the existence of two thresholds, one for the small-world transition and 
the other for the enhancement of synchronizability, in SW networks using a more rigorous analysis. They 
obtained that the SW regime starts at Si ^ 1/N whereas the threshold beyond which the synchronizability 
is improved goes like S'sync ~ k. This means that when k is large, Aj^^ of the underlying regular network 
contributes significantly to the synchronizability, and then it can be enhanced without additional shortcuts. 
This is manifested in Fig. [T^] by the fast decrease of S'sync when k approaches the critical value k^ . On 
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Figure 12: (color online) Synchronizability thresholds Ssync (o) for graphs with A'^ nodes (A'^ = 300, 400, 500, 1000) and 
range k S [1,70], averaged over 1000 realizations. Solid lines are based on an analytical perturbation (Eqs. II69I70I I valid in 
N^/^ < k < fcsync)- ^'^^ most parameters, Ssync lies within the SW region between the dashed lines (depicted for A'^ = 1000), 
but it is distinct from its onset Si. Inset: decay of the average distance £, clustering C, and eigenratio (squares) as shortcuts 
are added to a regular network of n = 500 and k = 20. We define Si and Sq as the points where £ and C are 75% of the 
regular network value. From [9l| . 
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Figure 13: Eigenratio R (left), distance £ (right) and eigenvalue A2 (inset) as functions of k (the range of nearest neighbours 
in the WS model of SW networks). N = 500 nodes, averaged over 100 realizations. The average number of shortcut per node 
S is fixed, so that the rewiring probability is p = S/2k; S = 0.3 (circle) and S = 0.5 (square). 



the other hand, at low k, the approximation in Eq. ([^^ is not vahd. However, the results in Fig. fT^ show 
that iSsync £ [0.3, 1] depends on k, but not noticeably on N. 

One observation from Fig. [T^] is that smaller distances are not be necessarily correlated to enhanced 
synchronizability as intutively believed. Indeed, if we keep the number of shortcut per node fixed, 5 ~ 1, 
and increase k as in Fig. 1121 we can see that the network distance decreases monotonically, while the 
eigenratio does not always follow the same trend. There is a range of k where the synchronizability of 
Type I is reduced (R increases) while the network distance becomes smaller, see Fig. [131 The Type II 
synchronizability (A2, unbounded MSF), however, is enhanced. In summary, it is not very meaningful to 
compare the synchronizability of two SW networks (with different N, k or p) considering only £. 

The relationship between £ and the synchronizability of the system was also scrutinized for SF networks 
in [l03l | , an important work that raised the interest of studies on the structure-synchronizability relation in 
complex networks. As seen in Fig. 1141 for random SF networks, i decreases when the degree distribution 
becomes more heterogeneous (decreasing of the exponent 7), however, the network becomes less synchroniz- 
able, since R increases. In these simulations, the mean degree of the network (fc) = fcmin^E^ also changes 
with 7, as well as the standard deviation of the degree distribution. As will be clarified later on (e.g., 
Eq. in this case, the eigenratio R is controlled by the heterogeneity of the degree distribution, being 

R ~ ^max/^min if ^min IS lar ge enoug h, w hile t he eigenvalue A2 has a dependence on the mean degree when 
the minimal degree is fixed 111 , 112{ . In [l03j the authors observed similar results in a SW network model 
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Figure 14: Synchronizability of SF networks of size = 1024. The average network distance (a) and the eigenratio (b) for the 
random SF model with fcmin = 5. The inset of (a) shows the mean degree (fc) and the standard deviation s of the connectivity 
distribution. Inset of (b): the maximum normalized load bmax. The horizontal lines in (a) and (b) indicate the values of £ and 
the eigenratio R for the 7 = oo case. The solid curve in (b) is the lower bound fcmax/femin in Eq. H79| l. The upper bounds of 
in Eq. I I79I I arc above the limits, but follow the same trend. All quantities arc averaged over 100 realizations. From |103|| . 



with hubs, and they considered them counterintuitive. Similar to the SW networks with different range k 
in Fig. [T31 where the eigenratio R is not simply controlled by £ , here it is also not surprising to observe that 
synchronizability can be reduced when £ becomes smaller [103|. 

Besides, it should be possible to observe the situation that R decreases at smaller i in random SF network 
for suitable combinations of parameters N, 7 and fcmin, for example, when fcmin increases at fixed 7 and N. 
Therefore, the conclusion is that the synchronizability of complex networks cannot be assessed solely on the 
average shortest path length £. 

Synchronizability dependence on betweenness centrality. In fl03!| it is argued heuristically that this appar- 
ently surprising behavior (smaller £, less synchronizability in SF networks) is due to the fact that a few 
central oscillators interacting with a large number of other oscillators tend to become overloaded. When 
too many independent signals with different phases and frequencies are going through a node at the same 
time, they can cancel out each other, resulting in no effective communication between oscillators. Thus 
the authors were motivated to examine the influence of the load (betweenness) of the nodes. It was shown 
(see the inset of Fig. [H] (b))that the synchronizability follows the same trend as the maximum load ^max 
(normalized b y the tota l load of the network). However, the load of a node in SF networks is closely related 
to the degree |Ti3l.lTil. i.e., nodes with large degrees or links connecting nodes with large degrees have, on 



average, a large load. The correlation between reduced synchronizability and heterogeneous load has also 
been observed in a variant of the WS network mo del b y adding m shortcuts to the network from a randomly 
selected node to one out of the Hc center nodes |l03| . In this case, when Hc is small, the m shortcuts are 
connected to a few hubs, and the degree becomes more heterogeneous and the maximum load &max increases, 
while the synchronizability decreases. As before, in these two examples, it is not very clear whether the 
change of synchronizability is mainly influenced by the degree heterogeneity or by the load itself, because 
these two properties are closely related. 

In the original WS network Q , the maximum load 6max decreases when the degree distribution becomes 
more heterogenous as p is increased. Based on this observ ation , it was claimed that more homogeneous 
load predicts better synchronizability on complex networks |lOCll |. Howe ver, a direct relationship between 



load and synchronizability can not be clearly established. In fact, in 115l | it is shown an example where the 
network displays improved synchronizability while the load becomes more heterogeneous when an original 
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Figure 15: (color online) The relationship between graph distance i and clustering coefficient C. Left: the original networks 
are the SW networks. Right: the original networks are the extensional BA networks. The different symbols are for different 
heterogeneity of degrees, meas ured by the standard deviation a of the degree distribution. All the data are the average over 
20 different realizations. From [llSj . 



random SF is rewired to obtai n non trivial clustering. 

The heuristic argument in |l03t is not clearly justified. If the picture described is correct, i.e., signals 
can cancel out each other at the central nodes, resulting in no efi^ective communication, then the Type II 
synchronizability (A2, unbounded MSFs) should also be reduced significantly. However, as seen in Fig. [TTl 
A2 in SF networks is mainly determined by the minimal degree fcinin, without a noticeable dependence on the 
degree distribution exponent or load. Moreover, when the minimal degrees becomes larger at fixed 7, the 
maximal load increases, however, both types of synchronizability are enhanced (smaller R and larger A2), 
in contrast to the heuristics. In fact, several further investigations have shown that the highly connected 
oscillators synchronize faster among the m an d form synchronization clusters [48j, Ill6l . IIITI . [SOj , which is in 
contrast with the argument provided in |l03| . 

Synchronizability dependence on the clustering coefficient. In [s^ it is pointed out that , in general, the 
eigenratio R increases with increasing clustering in a modified version of the BA model 118l |. Unlike the 



original one, in this model, motivated by the evolution of language, a new node is first linked to an existing 
node according to the preferential attachment rule, and also linked to the neighbors of this target node. 
Thus this model displays nontrivial clustering while keeping the same SF degree distribution. Besides, the 
eigenratio R is larger than that of the BA model. 

A recent work demonstrated that for both SW and SF networks, large values of clustering hinders global 
synchronization of phas e osci llators, since the network splits into dyna mica l clusters that oscillate at different 
frequencies [41, 43|. In [llSj . using the rewiring scheme proposed in |ll9| . that changes the clustering but 
keeps the degree sequence, it was shown that the eigenratio R also increases with C, for both SW and SF 
networks. This indicates that synchronizability is reduced when clustering C increases. 

Note that in the above investigations different structural properties change at the same time 91, lOd 

Il03 l| when the parameters characterizing the original networks are modified. Once again, it is difficult to 
draw conclusions about the relationship between one single st atisti cal descriptor of the network and its 
sy nchro nizability. Special attention was paid to this problem in 115l | , showing that in the rewiring scheme 
of [ll9l |. £ is correlated with the clustering coefficient C (Fig. [T5|). 

Synchronizability dependence on degree correlations. In many real-world networks, the degree of a node is 
often correlated with the degree of the neighboring nodes. Correlated networks show assortativ e fd isas- 
sortative) mixing when high degree nodes are mostly attached to nodes with high (low) degree [l20|. In 
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Figure 16: (color online) Synchronizability of degree correlated SF networks of size TV = 1000, fcmin = 5 and 7 = 2.5. Left: 
Behavior of th e eig enratio R. Right: behavior of the second lowest eigenvalue A2. Both as functions of the correlation coefficient 
rfe defined in '120||, for 7 va ryin g from 2 (blue line) to 5 (red line) in steps of 0.2. The results have been averaged over 100 
different realizations. From 11221 . 



practice, the degree-degree correlation of a network can be calculated as the Pearson correlation coefficient 
between degrees {ji, ki) of the nodes at the ends of the Ith link, i.e., 



ink) ~ ih) 



2 



- {kf) - iky ' 

where (•) denotes average over the total number of links in the network. Positive values of indicate 
assortative mixing, while negative values refers to disassortative networks. A spe cific correlation can be 



obtained by rewiring the links while keeping the degree sequences unchanged 121|. However, one can expect 
that i and 6max change also when the networks are rewired. 

The influence of degree correlations on synchronizability was addressed in js^] , showing that the eigenratio 
R increases when increasin g the assortativity of the network. A more detailed and systematic analysis 
was later on carried out in [l22| . They generated random SF network as in [l02| with a minimal degree 
fcmin = 5 and obtained the desired degree correlation using the rewiring procedure of Refs. jl2ll . l33 |. The 
dependence of the synchronizability on rk is shown in Fig. 1161 The main effect on the eigenratio R comes 
from the fact that A2 decreases when rk grows, while Xn remains roughly constant |l22l |. As it happens for 
the other parameters, the dependence of A2 on rk can be obtained from graph theoretical analysis, which is 
the subject we are going to discuss next. 

4- 1.5. Graph theoretical bounds to synchronizability 

Previously we have seen that many structural properties can influence the synchronizability of the net- 
works, but none of them can be regarded as the exclusive factor in the observed dependencies. The moral 
is that all works described in the preceding paragraphs seem not to be on the most appropriate way to 
elucidate the dependence of synchronizability on the network characteristics. Since the synchronizability 
depends on the two extreme eigenvalues of the Laplacian matrix, a sound analysis must attack the raw 
problem of the spectral properties of networks from a mathematical point of view, given that the simulation 
experiments are far from being conclusive. 

Graph theoretical analyses of the Laplacian matrix L, in the context of synchronizability, mainly focus 
on the bounds of its extreme eigenvalues. The implicat i ons of these bounds on differ ent types of complex 
networks have been discussed in several works 123 , 106 , 124l Il25l . Il26l . llOSl . Il04l . Il27l | . Here we summarize 



the main results and discuss how they help to understand the synchronizability of complex networks. Some 
detailed pro of o f the results c an be found in the corresponding references and in monographs on graph 
theory, e.g. [Hi IHi, EM \lE^ ■ 



First, we discuss the bounds for networks with pre scribed de gree sequences fcmin = fci<fc2<---<fcAf = 
fcmax- The bounds satisfy the following relations (see 128 . 12^ ) 



2 (1 - cos(^)) e(G) < A2 < ^^fcmin, (74) 

and 

TV , ^ 

N < 2fcmax, (75) 
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where e(G) is the edge connectivity of the graph, i.e. the minimal number of edges whose removal would 
result in losing connectivity of the graph. It follows that 

k~- - (l-cos(^))e(G)' ^ ' 

Note that the edge connectivity e{G) < fcmin- The equality holds only for special cases, in particular for 
networks which are homogeneous in degree. Then for a network with a fixed number of nodes and edges the 
edge connectivity ranges from 1 (a structured network with communities connected through single links) to 
a maximum value for a homogeneous network, and this turns out in the large variability of the lower bound 
for A2. 

The upper bound in Eq. (|74p is approached when the network is random. As shown in lOd . 124j . in a 



fc-regular random network, where each node is randomly connected to other k — fcmin nodes in the network, 
X2 = k — Oi^/k) as N—^_oo. In fact, A2 = Ckk, where Cfe — > 1 as k — > 00 both for /c-regular networks and 



random SF networks [106[. In this sense, the observation in jl05l | that A2 is almost constant, practically 
unrelated to fcmin, is incorrect. So for large random networks with large enough minimal degree, A2 ~ fcn 
is a good approximation. 



Another general lower bound for A2 is [132 1 



N 

< A2 < -fcmin, (77) 



so that 



ND - N-1 



<i?<^^-^^, (78) 



k ■ ~ 

'''mm 



where D is the diameter of the graph, i.e. the maximum value of the shortest path lengths between any two 
nodes. 

Let us discuss the implications of the above bounds for WS networks. In this model, when the probability 
of having shortcuts is very low, -pkN <C 1, then D ^ N /k and the lower bound in Eq. ([77|) approaches zero 
as (regular network). Beyond the onset of the SW regime {pkN ^ \), D decreases and approaches 

D ^ In(iV), and A2 increases for fixed N. Thus this bound allows us to understand why the network 
synchronizability is inversely proportional to as observed numerically in [oilflOQi fFigs.[51and [TU)) . However, 
A2 is not immediately bounded away from zero, since it approaches zero faster than 1/N just after the onset 
of the SW regime. When moving deeper into the SW region, pkN ^ N so that each node has S ^ \ 
shortcuts, and A2 is already bounded away from the lower bound A/ND and approaching the upper bound 
fcmin- In this regime, A2 will not be sensitive to changes in the diameter D. This helps to understand why 
the synchronization threshold is different from the onset of the SW behavior in Fig. [T% 

Thus, in WS networks with high p, fcmax/fcmin provides a good lower bound to the eigenratio R. The 
upper bound, iVZ?fcmax/2, still increases with N, and can be several orders of magnitude larger than the 
actual value of R even for small random networks. In fact, R may not follow the variation of the distant 
upper bound when the diameter D or the size N change. Thus, in complex networks with both local and 
rando m co nnections, a close relationship between the synchronizability and £ is not expected. 

In I103I similar bounds are obtained but as a function of £ and &max as 



^ <R< Nk^^^bm^^De. (79) 

"'min 



These authors pointed out that experimental values of R are closer to the lower bound, and far away from the 
upper bound (Fig. I14p . Thus quite probably such upper bound does not provide meaningful understanding 
of the relationship between synchronizability and £, D or 6max, since the change of the upper bound by 
these structural measures is likely not to be reflected in the actual value of R. Indeed, for the SF model 
considered in this paper and arbitrary random enough networks with suffiently large mimimal degrees, recent 
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Figure 17: A schematic plot of the bounds of A2 for networks with minimal degree k„ 



results pjjj showed that the eigenvalues are bounded as 

2 

fcmin(l 7=) < A2 < k„ 

V W 



2 

fcmax < Aa, < /Cn,ax(l + -7=), (80) 

V w 

which can explain more clearly how the eigenvalues in Fig. Illla nd in ^lOS] depend precisely on the maximal 
degree fcmax and the mean degree (fc) . A recent advance [112| in spectral analysis of random SF networks 
shows that the ma ximal eigenvalue is very close to the lower bound in Eq. ([75|) . Ajv — fcmax + 1- The 
observation in |l03l | that the synchronizability (Type I) is anti-correlated with a more heterogeneous load 
distribution, and smaller distance, is a consequence of the positive correlations between the load and the 
degree, and negative correlation between the distance and the degree, in the models considered. Indeed, 
when keeping the degree sequence unchanged, it is observ ed th at more heterogeneous load distributions can 
be positively correlated with synchronizability of Type I [llSj . 

The above analysis of bounds provides justification about why we can observe different R for increased 
or decreased heterogeneity, distances or loads. Moreover, the bounds expressed by these quantities are not 
tight at all in the particular examples, e.g., in Eq. ([7^ . In general, the graph theoretical analysis states 
that randomness improves the synchronizability, since A2 is well bounded away from 0, while in networks 
with dominantly local connections, A2 approaches to in large networks. In other words, for a prescribed 
degree sequence, the eigenratio R changes mainly because of A2. A schematic plot of the bounds of A2 is 
shown in Fig. [171 

On the other hand, the local organization of the connections even within a small part of the network, 
can make the eigenvalue A2 deviate significantly from the upper bound ki-ninN/{N — 1). For an arbitrary 
network, there is a better upper bound for A2 as 

A2 < 2ig, (81) 
where ig is the isoperimetric number of a graph [l33l Il25l . Il26| , which is defined as 

.,=min^. (82) 

Here S C G is a, subset of the nodes, with Q — S denoting its complement, and |5| is the number of nodes 
within S, with < |5| < N/2. Besides, \dS\ is the number of connections between S and its complement, 
namely, 

\9s\ = Yl E (83) 



For an arbitrary partition of the network into S and Q, we have [12; 



(84) 



This means that the synchronizabihty of the network is determined by the sparse connections between the 
two subnetworks. For instance, if a smaU set S made up of, say, 20 nodes, is connected to Q ~S with just one 
Unk, then A2 < 0.1, regardless of how the nodes are connected within S and within the large complement 
Q — S. It follows that the statistical properties of the network Q are mainly determined by the huge part 
G — S, while A2 is independently constrained by the small subgraph. This result is in complete agreement 
with the path towards synchronization in modular networks presented in scction [3.1.51 where the community 
structure has been demonstrated to impose scales in the synchronization process. 

Everything up to now indicates that statistical properties, such as degree distribution, i, etc., may 
not be always correlated with the synchronizabihty of the network. In fact, it was particularly shown in 
124Lll25l . ll26l |. that networks with very different synchronizabihty can be constructed for the same prescribed 



degree sequences, because ig can be at any place in a broad range between the lower and upper bounds in 
Eq. JTll), see also Fig.fTTl 

Another similar bound in graph theory is based on the Cheeger inequality 134L 122 |. 



X2<hg, (85) 

where hg — min^ hg{S) and 

- MI- <'^' 

As for correlated networks, in [l22] the authors showed that in random networks with degree correlations 

dhg{S) 



dr, 



k 



< 0, (87) 



which explains why the synchronizabihty is reduced in networks with assortative connections, as shown 
in Fig. 1161 Along the same line, the dependence of the minimum nonzero eigenvalue on the topo logical 



properties of the network and its degree-degree correlation coefficient r was also analyzed in [127[. The 
authors derived a rigorous upper bound for A2 as, 

(fc)(fc^)~(fc3)(fc2) 
^^-(^-"^ (fc)((fc2)_(fc)2) ■ (88) 



It is worth mentioning also the inequality presented in [1331 . 112: 



A2 > A:„iax - y A:2^ax - ig- (89) 

Note that k^^x — \J ^max ~ *g 9- decreasing function of fcmax if ig is fixed. Based on this bound, we can 
expect some apparently counterintuitive effects on the synchronizabihty of complex networks. Suppose we 
put more links to the graph t/, but only add them to the nodes within Q — S and S, but not to the nodes 
between S and Q — S so that ig does not change. For simplicity, we can further assume that the nodes 
within S and Q — S are connected with a uniform probability / (random networks), so that fcmax and the 
mean degree (fc) increase when more and more links are added. In this case, the synchronizabihty of the 
subnetworks S and Q — S is enhanced at larger /, see Fig.[9l However, A2 of the whole network Q is reduced 
according to Eq. (|89p . Therefore, in the case of two coupled networks, enhancing the synchronizabihty of 
the subnetworks may actually reduce the synchronizabihty of the whole network. Phenomenologically, this 
is intuitively expected, because the subnetworks tend to form distinct synchronized clusters. 

Based on the above arguments (Eqs. ()81l85l89p ). networks possessing a clear community organization 
display a small synchronizabihty, since the density of conne ction between different communities can be much 



smaller than the density within the communities |135l . [13; 



Recapitulating, we have seen that for a prescribed degree sequence, it is possible to construct a very 



large number of networks ranging from fully local connections to fully random networks |124| , with many 
possible structures in between. However, the degree sequence by itself is not sufficient to determine the 
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synchronizability. On the other hand, we have seen that the synchronizabiHty is not directly related to graph 
measures, such as distance, clustering or maximal betweenness. Admittedly, the weak connections between 
two subnetworks (characterized by the isoperimetric number ig) determine the behavior of the eigenvalue 
A2, and hence that of the synchronizability of the whole network. The bounds discussed in this section are 
valid for any network. However, we would like to point out that one needs to be careful in the interpretation 
of these general analytical results. A linear relationship between the bound and some network descriptor 
does not mean that we can always expect the same relation to hold between these descriptors and the actual 
eigenvalues for particular types of networks. The reason is that the bounds, the network descriptors and 
the eigenvalues of the Laplacian can follow totally different scaling laws in particular types of networks. 
For example, for random SF networks with large enough minimal degree, R ~ fcmax/fcmin ~ A^^/(^~^) from 
Eg . [501 while the upper bound in Eqs. (|78l79p increase faster than N ■ A^i/(t~i). The numerical results 
in [103| showed that the eigenratio R and the load follow different scaling laws also. This indicates that 
one must be cautious not to generalize the observation of such correlations in one type of network, not to 
interpret such observed correlations as the ultimate responsible for the synchronizability without a deep 
analysis of their constrains. 



4-1.6. Synchronizability of weighted networks 

Up to now, we have considered the influence of the network topology on synchronization, assuming that 
the connection weights are the same for all the links in the network, i.e., the networks arc unweighted. 
However, this is not the case for many real- world networks. Indeed, many complex networks where synchro- 
nization is r elevant ar e actually weighted and display a highly heterogen eous distribution of both degrees 
and weig hts [133, im, [isi, Il40| | . Examples include neural networks [lij, Il42| , ai rport networks (137.] and 



the structure of the networks characterizing epidemic outbreaks in different cities 143 . 144j . Furthermore, 



it has been observed that in many cases, the connection strength is not an independent parameter, but it 
is correlated to the network topology. The analysis of some real networks [137] yields the following main 
properties: 

(i) the weight Wij of a connection between nodes j and i is strongly correlated with the product of the 
corresponding degrees as (wij) ~ {kikj)^; 

(ii) the average intensity S{k) of nodes with degree k increases as S{k) ^ k^ . Here the intensity Si of a 
node i is defined as the total input weight of the node: 



N 

J^azjW^j. (90) 



Note that the inclusion of a distribution of weights in the network affects directly its classification within 
topological homogeneity or heterogeneity. For example, a regular lattice with a very skewed distribution of 
weights can eventually represent a SF topology. From a mathematical point of view, the adjacency matrix 
is in this case simply substituted by the weight matrix. On the contrary, from a physical perspective, it is 
still interesting to keep separated the topology of interactions from the distribution of weights, and answer 
questions, whenever possible, discriminating these two topological aspects. 

The first works on synchronization in weighted networks considered that the weighted input of a node i 
from a node j depends on the degree ki of node i 9^ 3^], with a model of weighted coupling as Wij — kf, 



so that the matrix G — (Gij) in Eq. ((55)) reads 

G^D%D-A) = D''L. (91) 

Here Dij = Sijki is the diagonal matrix of degrees. is a tunable parameter that keeps the network topology 
unchanged, but varies the distribution of the weights of the links. 

Within this scheme, the weights between a pair of nodes i and j are in general asymmetric, because 
Wij = kf and Wji = kj. However, since 

detiD^L - XI) = detiD'^/^LD"^^ - XI), (92) 
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Figure 18: Eigenratio R as function of 9 in Eq. I|91ll : (a) random SF networks with 7 = 3 (•), 7 = 5 (H) and 7 = oo (solid 
line), for fc^jn = 10; (b) random networks with expected SF sequence for 7 = 3 and fc^in = 10; (c) growing SF networks for 
7 = 3 and m = 10; (d) SW networks with M = 256 (•) and M = 512 (■), for k = 1. Each curve is the result of an average 
over 50 realizations for N = 1024. Modified after |32|| . 



is valid for any A, the spectrum of eigenvalues of the matrix G is equal to the spectrum of a symmetric 
matrix defined as if = D^^^LD^^^. As a result, all the eigenvalues of G are real, and the synchronizability 
can still be characterized in the framework of the MSF. 

The synchronizability of various complex networks as a function of the parameter 9 is shown in Fig. 1181 
Except for k-regular networks, in all other cases, including the SW networks, the eigenratio R exhibits a 
pronounced minimum at 9 = —1. Here the SW networks are obtained by adding M < N{N — 2k — l)/2 
new links between randomly chosen pairs of nodes on the basic regular array where each node is connected 
to its 2k first neighbors. 



In [95|, I34I the authors also characterize the synchronizability of the network, related to A2 , using the 
notion of the cost of the network. When Eq. (|64p is satisfied, the fully synchronized state is linearly stable 
for a > CTinin = The cost is defined as the total input strength of the connections of all nodes at the 

synchronizability threshold: CTmin j Wijttij — (Tmin X]i=i ^i- ^ more convenient definition for comparisons 
is obtained normalizing by the number of nodes, such that 

Co ^ = (5)/A2, (93) 

IN ai 

where {S) is the average intensity of nodes in the network. Similar to R, Cq is also minimal a,i 9 — —1 

(Fig.IHl). . rnrrh . 

Interestingly enough, in [95, 32"! it was obtained that in SF networks with fixed minimal degrees fcmin, the 
weighted versions {9 = —1) behave differently to the unweighted networks when one looks at the dependence 
of both the eigenratio R and the cost Co on the scaling exponent 7, as shown in Fig. [2D] . 

= —1 is a special case. The coupling matrix is now G = D~^L, and all the diagonal elements Gu = 1 . 
It is usually called the normalized Laplacian of a graph. Based on graph spectral analysis results in \l2^ 
for random networks with arbitrary given degrees, it can be shown that the spectrum of the normalized 
Laplacian tends to the semicircle law for large networks. In particular, for fcniin 2> -y/ (k) N, one has 



max{l-A2,Ajv-l} = [l + 0(l)]^=. (94) 



This result is rigorous for ensembles of networks with a given expected degree sequence and sufficiently 
large minimum degree femin, but the numerical results reported in Fig. [20l support the hypothesis that the 
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Figure 19: Normalized cost Co as a function of in Eq. H91|l for random SF networlis with 7 = 3, fcmin = 10 and = 1024. 
Modified after 
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Figure 20: Eigenratio _R (a) and normalized cost Co (b) as a function of 7 for random SF networks with with 6 = (o) and 
8 = —1 (•) and for random homogeneous networks with the same average degree (fc) (o). The other parameters are fcmin = 10 
and N = 1024. The dashed line corresponds to 7 = 00 ((fc) = 10). The solid line in (a) is the approximation given by Eg. 1951 
Modified from 13211. 
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approximate relations 

A2 w 1-2/7(1)", Aat w 1 + 2/7(1)" (95) 

hold under much milder conditions. In particular, the relations ()95p are expected to hold for any large 
network with a sufficient number of random connections, fcmin 3> 1. 

Furthermore, the synchronizability in this case seems to be independent of the degree distribution. It is 
only controlled by the average degree (fc) , since the synchronizability of the weighted SF networks is almost 
identical to that of a regular random network where each node has the same degree ki = (k) . These results 
demonstrate that the topological degree of networks is not the only determinant of the synchronizability of 
the networks; having a heterogeneous distribution for the connection strengths can significantly influence 
the synchronizability. 

4-1.7. Universal parameters controlling the synchronizability 

What are the leading parameters governing the synchronizability for the more general case in which 
weighted networks are considered? It [111.] it has been proved analytically and verified numerically what 
controls the synchronizability of sufficiently random networks with large enough minimal degree (fcmin S> 1). 
It is the distribution of the intensities Si defined in Eq. ([50]) . The intensity of a node incorporates both the 
information about the topology and the weights of the connections in the networks. The main finding is 
that the synchronizability is sensitively controlled by the heterogeneity of the intensity Si. The eigenratio 
R and the normalized cost Cq can be expressed as 

R = An^^ Co^Ac^, (96) 

where S'min, S'max, and (S) are the minimum, maximum and average intensities, respectively. The pre-factors 
Afi and Ac are expected to approach 1 for large average degree (fc) . Equations are universal in the sense 
that they apply to many random networks with arbitrary degree and weight distributions provided that the 
minimal degree is sufficiently large. The main hypothesis behind this result is the assumption that the local 
mean fields Hi = (l/fcj) J2j '^jjH(xj) can be substituted by the global mean field Hi w H = (l/N) H(xj). 

In Fig. Eq. is corroborated by numerical results of R and Co for networks with several degree 
and intensity distributions of degrees, using the weighted coupling scheme 

= Si/ki. (97) 

This coupling scheme means that the intensity of the nodes, which is not necessarily correlated with the 
degrees, is uniformly distributed into the input links of the nodes. It covers the coupling scheme in Eq. I|9ip 
as a special example: Si = k]^^ . Recent analysis in [l45j on the spectral density of SF networks with a 
weight ed L aplacian matrix similar to Eq. (^1]) . also confirms that Eq. ((5^ holds. 
In lllj the authors also presented results for the following coupling scheme 

w^j = {hkj)\ (98) 



that describes the relationship between the weights and the degrees in some real networks [1371 Il46l | . The 



tunable parameter 9 controls the heterogeneity of the intensity Si and the correlation between 5*^ and fci, 
since Si = kl'^^{kj)i, where (fcj')i = (l/fci)X]^| is approximately constant for ki ^ 1 when the degree 
correlations can be neglected. Variations of 9 have a significant impact on the synchronizability of networks 
which are heterogeneous in degree(Fig. [25]). Note that in heterogeneous in degree networks, the weighted 
coupling in Eq. (|98p may result in a broad distribution of the input weights Wij among the ki links of the 
node i, especially when 9 is not close to 0. However, as shown in the insets of Fig. [221 for various networks 
and 9 values, R and Co collapse again to the universal curves when regarded as functions of Smax/S and 
(S) /Smin, respectively. The fact that the universal formula holds for a broad range of 9 values shows that 
the mean field approximation used to obtain Eq. (|96p often remains valid under milder conditions. 

It is important to stress that these results also hold for unweighted random networks. In this case 
Si = ki for all the nodes and one gets the results in Eq. ([80]) for large random networks with minimal degree 
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Figure 21: (a) R as a. function of 5max/Sniin and (b) Co as a function of {S)/Siain, averaged over 50 realizations of the 
networks. Filled symbols: uniform distribution of Si G [Smim S'lnax]- Open symbols: power-law distribution of Si, P{S) ^ 
for 2.5 < r < 10. Different symbols are for networks with different topologies: BA networks (o), growing SF network with 
aging exponent a = —3 (□), random SF network with 7 = 3 {o), and fe-regular random networks (A). The number of nodes 
is N = 2^^ and the average degree is (fc) = 20. Insets of (a) and (b): Aji and Aq as functions of (k) for Smax/Smin = 1 (°), 
2 (□), 10 (A), and 100 (*), obtained with uniform distribution of S j in fc-regular networks. The dashed lines are the bounds. 
SoUd Unes in (a) and (b): Eqs. ^ with Ar = Ac = 1. From [uj]. 
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Figure 22: (a) Eigenratio R and (b) cost Co as functions of 9 for BA networks (o), growing SF networks with aging exponent 
a = —3 (□) and random SF networks with 7 = 3 (A). Each symbol is an average over 50 realizations of the networks with 
(fc) = 20 and A'^ = 2^*^. Inset of (a): the same data for R as a function of S max /Smin ■ Inset of (b): the same data for Co as a 
function of (S'>/5niin- Solid lines: Eqs. (|96)l with Ar = Ac = 1. From p^ . 
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Figure 23: (a) Eigenratio R as a function of Smax/S'niin for BA growing networks with m = l((fc) = 2) (o), m = 2{{k) = 4) 
(□) and m = 10{{k) = 20) (A), (b) R for SW networks with (k) = 20 and rewiring probability p = 0.01 (o), p = 0.1 (□) and 
p = 0.3 (A). SoUd lines: Eqs. (|96]l with Ar = Ac = 1. 



^min ^ 1- These resuhs provide much tighter bounds than those discussed in the previous sections (e.g., cf. 
Eqs. (|76I78[[7^ ) which depend on the system size TV. Interestingly, Eqs. ([M)) also provide meaningful insights 
into the problem for other special networks. For example, consider the class of SF networks generated using 
the BA model. When ni — 1, the network is a tree, and R is much larger than S'max/5'miii- However, it 
is an increasing function of S'max/'S'min, showing that the heterogenity of the intensity is still an important 
parameter. But for m = 2 ((fc) = 4) it approaches the universal curve quickly (Fig. [21] (a)). The drastic 
change of synchronizability from m = 1 to m = 2 can be attributed to the appearance of loops [89] . In WS 
networks [9| with N = 2^° and (k) = 20, R collapses to the universal curve even when the networks are 
dominated by local connections, e.g. for a rewiring probability p = 0.3 with intensities S'max/'S'min ^ 10, see 
Fig. [231(b). 

4--2. Design of synchronizable networks 

An interesting subject related to the impact of network structure on synchronization dynamics is the 
design of synchronizable networks. Here we review several ideas exploring this issue: weighting the couplings 
leaving the topology unchanged, perturbing part of the network topology, and finally searching for optimal 
topologies with respect to synchronizability. Note that the following theoretical schemes may not directly 
apply to real complex networks. It is difficult to conceive real systems where the weights can be tuned at 
discretion, or where the topological substrate of interactions can be changed accordingly. Nevertheless, the 
insights given by these works allow for a deeper understanding of the synchronizability of networks. 



4-2.1. Weighted couplings for enhancing synchronizability 

The previous analysis shows that for networks that are heterogeneous in degree, synchronizability can 
be enhanced by balancing the heterogeneity in the degree distribution with suitable weighted couplings, 
towards the obtention of a homogeneo us di stribution of the intensities 5"^ 

A different scheme is presented in 



147l |. where it is assumed that the weight of a link is related to its 



betweenness ba as 



b"- 



E 



6" 



(99) 



where a is a tunable parameter that controls the dependence of the weights Wij on the loads 6y . The 
zero-sum requirement of the matrix G implies that Ga = 1 for all i. Note that a = corresponds to the 
weighted coupling scheme in Eq. ([OT|) at the optimal point 6 = —1. As seen in Fig. [231 the eigenratio R 
depends on a, reaching a minimum at a value < a < 1, showing that the synchronizability in SF networks 
can be slightly enhanced compared to the optimal case of the weighted coupling scheme [H, llllj . 
The SF netw ork c onsidered in 147l | is a generalized BA model with a preferential attachment probability 



TTi ^ ki + B 
distribution, 



148 1 



and ' 



where the parameter B controls the exponent 7 = 3 + B/m of the power law degree 
'' = fcmin- At large a values, only the links with the largest load s are significant, which 



can lead to effectively disconnected nodes, so that synchronizability is reduced. In jl47l | it is also pointed out 
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Figure 24: (color online) (a) eigenratio \m/^2 (in logarithmic scale) for SF networks in the parameter space {a,B). (b)The 
relative synchronizability F = log{Xpf / X2) — [log(A]v/A2)]a=o vs {a,B). In all cases m = 2, and the graphs refer to averaging 
over 10 realiz ation s of networks with N = 1000. The domain with with F < is outlined by the black contours drawn on the 
figure. From |147|| . 



that for large minimal degrees, the regimes corresponding to enhanced synchronizability are reduced so that 
the minimum approaches to a = 0. This demonstrates again that for random networks with large enough 
minimal degree, Eq. (P5|) is asymptotically valid regardless of the detailed weighted scheme, as claimed in 



ml. 



The explanation of the observed enhanced synchronizability proposed in is that the load bij reflects 
the global information of the network, while at a = only the local information (degree) is employed. 
Such a heuristic explanation, however, is not supported by several further investigations. In fact, only the 
local information can also lead to similar enhanced synchronizability. For example, consider that the weights 
depe nd on the degrees following Eq. ()98p , and then normalize to allow fully uniform intensity Si = I, namely 

m 



G. = - ^' '7, ... =-^^^- (100) 



Aga in, a = co rresp onds to the optimal case (Eq. ([9T|) . 6* = — 1) of the weighted couphng scheme 9^ H, 



111] . Similar to 147|, the synchronizability can be further enhanced in a range a > (Fig. [55]). However, 
the minimum moves closer to a = when the networks are larger, indicating that the synchronizability 
in large random networks with femin 3> 1 can hardly be enhanced further with other weighted coupling 
schemes different from the optimal case in Eq. This coupling scheme was also considered in [t^ where 

synchronizability in SF networks of maps is enhanced for a > 0. Additionally, we n ote t hat the coupling 



form in Eq. (jlOOp has been recently revisited from the viewpoint of gradient-network [l& 
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Figure 25: Eigenratio R (a) and n orm alized cost Co (b) as a function of the parameter a in Eq. I|100| l. for random SF networks 
with 7 = 3 and fcmin = 10. From Il49ll 



Another work 15ll | introduced an additional parameter f3 into the couphng scheme in Eq. (jlOOp . as 



In this case, the intensity of the node is Si — (X^jer ^f)^ ^ ■ P ^ 1; intensity is not uniform and 
becomes more heterogeneous when |1 — /3| increases. For any given a, the synchronizability is optimal at 
(3 = 1 where the intensity is fully uniform. Besides, for a fixed /3, there is also a value of a for which the 
best synchronizability is achieved. 

More weighted coupling methods have been proposed. In [97| the authors use the information about the 
age of the nodes in growing networks and introduce asymmetrical coupling between old and young nodes 
(first and latest nodes to join the network, respectively). Old nodes in general have large degree and young 
nodes small degree. The authors propose a connectivity matrix 

G„- = „ «»^-Qu (102) 

where 9^ = (1 — 9)/2 if the connection is from old to young nodes {i > j) and 8^^ = (1 + 0)/2 if the 
connection is from young to old nodes (i < j). The limit 9 = —1 {9 = 1) gives a unidirectional coupling 
where the old (young) nodes drive the young (old) ones. It was shown that in SF networks, synchronizability 
is enhanced when couplings from older to younger nodes are dominant {9 < 0). When large heterogeneous 
degrees between the old and young nodes occur, this scheme is quite similar to those in Eqs. (|99llOOI10ip . 
However, in this case the eigenvalues are complex, and both, the ratios of the real and imaginary parts of 
the eigenvalues were employed simultaneously to characterize the synchronizability. 

In spite of these numerical observations, a clear understanding about why the synchronizability is further 
enhanced with the various weighted coupling schemes (Eqs. (|99llOOI10ip ) has not been obtained. The same 
scheme as in (jl02p . but without the normalization (G^ = —aijQij) was considered in 'l52!|, and again 
the synchronizability is enhanced for 9 < 0. In this case, the change of synchronizability is mainly due 
to the heterogeneity of the intensity distribution Si, which becomes more homogeneous for 6* < because 
the old nodes are hubs and have most of the connections to young nodes. In this case, the universality 
in Eq. ([M)) should apply when the minimal degree is large enough. Finally, it has also been shown that 
a random distribution of the weights of the connections of regular networks, with only nearest neighbors, 
can also enh ance synchronizability. This fact is related to the effective presence of short cuts in terms of 
weights |l53l |. 

As it has been already pointed out, the particular scientific course of action taken by proposing different 
weighting schemes to enhance synchronizability is unfinished, and probably an unfruitful quest given the 
many possibilities of inventing new weights. A more rigorous analysis of the eigenspectra of general graphs 
beyond the already obtained bounds is absolutely required to boost this line of research. 

On the other hand, the above weighted couplin g sch emes are static. In many real-world systems, the 



network structure evolves and changes with time. In [ll7j the authors proposed a scheme that can adaptively 
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Figure 26: The average input weight V{k) of nodes with degree fc as a function of k for the Rossler oscillators (o) and the 
food- web model (•) (a) and its dependence on various parameters: m (b), TV (c) and p (d). Results in (b-d) are averaged over 
10 realizations of the networks with random initial conditions. For clarity, only the results f or th e Rossler oscillators are shown 
in (b-d) and are logarithmically binned. The solid lines in (a-d) have a slope —0.48. From IllTll . 



tune the correlation between the degrees of the nodes and the weights of the hnks as in Eq. ([M]) . with 
9 « —0.5, so that synchronizabihty can be significantly enhanced compared to the unweighted counterpart. 
The adaptation scheme is based on the local synchronization between a node and its ki direct neighbors in 
the network. Each node tries to synchronize to its neighbors by increasing the connection strength among 
them. By doing this, the coupling strength of the node i with its neighbors increases uniformly trying to 
suppress the difference A,; with the mean activity of its neighbors, namely, 

G,,{t)^a,,V,{t), V,=pA,/{l + A,), (103) 

where = |H(xi) — (l/ki) ayH(xj)|, and p > is the tuning parameter. Note that with this adaptation 
scheme, the input weight {wij = Vi) and the output weight {wji = Vj) of a node i are in general asymmetrical. 

This adaptive process was simulated using Rossler oscillators and a chaotic food web model on BA 
networks, and both the unbounded and bounded MSFs were considered. For the unbounded case, the 
system approaches complete synchronization when p > 0, while for the bounded case, one has that this 
happens for < p < pc, where pc depends on the particular oscillators, and on the system size N. In 
both situations, when synchronization is achieved, the adaptation process will lead to a weighted coupling 
structure where the input strength of the links of a node displays a power law dependence on the degree as 

V{k) - fc^ (104) 

with 9 ~ —0.5. The results for the unbounded MSF, which are rather robust to variations in network 
models, parameters and oscillator models, are shown in Fig. [26l 

The mechanism underlying such a s elf-organiz ation of the weighted structure is due to the degree- 
dependent synchronization difference A^ (ll6l lll7l| . Starting from random initial conditions on the chaotic 
attractors, both the local synchronization difference A^ ^ 1 and the input weights for each node increase 
rather homogeneously in the whole network, i.e., Wij = Vi{t) « jt. Now, the intensity of the node Si{t) = 
Vi{t)ki = pkit. Hence, nodes with large degrees are coupled stronger to the mean activity of their neighbors. 
As a consequence, after a short period of time the synchronization difference Aj for those highly connected 
nodes decreases, and the weights Vi of different nodes evolve at different rates and converge to different 
values. Once synchronization is achieved, the input strength Wij = Vi is small for nodes with large degrees. 

The adaptation process makes the intensity more homogeneous, so that it is expected that the synchro- 
nizabihty is enhanced. In Fig. [571 we show the eigenratio R, as a function of the network size N. There we 
compare the original unweighted network with two weighted networks after the adaptation. Suppose that the 
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Figure 27: The eigenratio R as a function of A'^ averaged over 20 realizations of the networks. The solid lines are power-law 
fitting. The weighted networks are obtained by the adaptive process with the conditions: M = 5, 7 = 0.002 with H(x) = {x, 0, 0) 
for the Rossler oscillators (□) and H(x) = {0,y,z) for the food-web model (A). The networks are synchronizable if R < Ra: 
Rossler oscillators, Ra = 40 (dashed line), food-web model, Ra = 29 (dashed-dotted line). From 117) 



largest network size synchronizable for the bounded MSF (i? = 02/01) is A^i for unweighted networks and N2 
for weighted networks obtained from adaptation. It then follow s tha t N2/N1 ~ (02/01)^/^^+^-' ^ (a2/ai)^ 
for power law degree distributions regardless of the exponent 7 |ll7l |. 



4-2.2. Topological modification for enhancing synchronizability 

Some authors have proposed to enhance synchronizability by perturbing the network topology. Based 
on the ar gument that heterogeneity in the betweenness distribution is related to poor synchronizability, in 



154l Il55l | it is proposed to modify the nodes or links with the highest maximal betweenness. As already 



noticed, in SF networks, the betweenness of a node and a link is strongly correlated with the degree fc^, 
and the product of degrees kikj of the two nodes at the ends, respectively. The perturbation proposed in 



154j consists of dividing the node with the highest degree into a group of several fully connected nodes 



and redistribute the ki links equally over the new nodes. Following this scheme, the synchronizability can 
be substantially enhanced by modifying a very small portion of the nodes. The enhanced synchronizability 
follows closely the reduced maximal degree in the networks. It was also shown that the average distance 
actually increases when the hubs are divided. In [l55l | the authors propose that the connections with the 
largest kikj are broken, and again the synchronizability can be enhanced by cutting a small fraction of 
links with high betweenness. These ideas can plausibly be implemented in technological networks, where 
the substitution of hubs by a core of nodes is possible. In this way, the redistribution of load will improve 
traffic, and as a by-product, the synchronizability. 



4.2.3. Optimization of synchronizability 

A more straightforward approach to the d esign is that of asking which are the best network architectures 
to get an optimal synchronizability. In |l5d | an optimization scheme (e.g. simulated annealing) combined 
with a network rewiring algorithm to minimize the eigenratio R is applied. In this case, the total number 
of nodes and links are preserved. The resulting networks, with optimal synchronizability, called entangled 
networks, are found to be very homogeneous in many topological measures, such as degrees, distance between 
nodes, betweenness etc. This result is quite relevant because it provides a null model that allows to compare 
the synchronizability of networks directly with its optimal counterpart. 

A similar optimization sche me w as applied to study the optimal synchronizability in networks with 
a preserved SF degree sequence jl57l|- In this case, the synchronizability can only be slightly enhanced. 
An interesting finding is that the optimized networks become disassortative and the clustering and the 
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Figure 28: (color online) Eigenvalue ratio R as a function of the number of algorithmic iterations. Starting from different 
initial con ditions, with A'^ = 50, and (k) = 4, the algorithm converges to entangled networks with very similar values of R. 
Prom [Tsl. 



maximal betweenness is reduced, which is consistent with the observed enhanced synchronizabihty obtained 
by changing these features. 

Thus, optimization schemes are helpful to identify meaningful topological features that correlate with the 
synchronizabihty of complex networks. However, such optimization schemes are computationally demanding 
when we deal with large networks. The development of analytical tools to attack this optimization problem 
is currently a major challenge in the subject. 

Furthermore, in these two studies, the underlying network topology is bi-directional. In many of the 
weighted coupling schemes previously mentioned, the conn ectio ns are effectively directed, since the coupling 
strength is asymmetrical for the input and output links. In [l58j the authors went to the extreme by imposing 
an unidirectional information flow so that the networks become optimally synchronizable with R = 1. For 
any topology, the maximally synchronizable network can be achieved by imposing that the network: (i) 
embeds a directed spanning tree, (ii) has no directed loop, and (iii) has normalized input strengths. With 
these conditions, the original networks are changed into feedforward networks without any feedback, and 
= Ai < A2 = • • • = Xn ~ A. Furthermore, the synchronization of the whole network is achieved in a 
hierarchical way. The conditions (i)-(iii) lead to the following path towards synchrony. There is a node in 
the directed spanning tree that has no input and acts as a master oscillator driving the network dynamics. If 
a = (tA, the oscillators that are just one hierarchical level below the master oscillator will synchronize. Then, 
the next lower level oscillators will also get synchronized and so on until the whole network reaches complete 
synchronization. Note that this happens for the entire range of the coupling strength where a = aX. 

Note that the optimization schemes discussed above have considered maximizing the Type I synchroniz- 
abihty by minimizing the eigenratio R under certain constrains. Such optimized networks, however, may 
not have enhanced synchronizabihty of Type II and smaller cost of synchronization that are associated to 
the eigenvalue A2. Furthermore, in these studies, the underlying network topology is un-directed. 

Is it possible to find the network structures that have the optimal synchronizabihty of both Type I and 
Type II and the smallest cost, among all possible network configurations? In 96, 158] an elegant answer is 
provided. They assume that MSF has a bounded convex region in the complex plane and denote ai and a2 
as the thresholds along the real axis. For any networks of oscillators, if synchronization manifold is stable 
for the coupling strength CTmin < f < CTmax, the synchronizabihty and cost of synchronization can be defined 
as 

n 

_ 0"max ^ 

*Jsyn — ; ^ syn '^min / ^ "^ij ■ 

f^min - . 

^,3 
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For real eigenvalues ordered as in Eq. 1591 we have 



^ syn — ■ T ; ^ syn — "7 / ^i- 

Cmin Ajv A2 

i— 1 

As in Eq. 1931 the cost is the minimal total coupling strength when the network is just able to synchronize. 
In [gij l it is proven that for any network 

Ct2 

Ssyn < —, Csyn > ai{N - 1). 

If the spectum of the coupling matrix satisfys 

= Ai < A2 = • •• = Aat = 1 (105) 
then the network will achieve the maximal synchronizability and the minimal cost, i.e., 

Ssyn — , Csyn — ai{N — 1). 

ai 

In 9^ 158l | such optimal networks are obtained by imposing an unidirectional information flow. For any 
topology, the maximally synchronizable network can be achieved by imposing that the network: (i) embeds 
a directed spanning tree, (ii) has no directed loop, and (iii) has normalized input strengths. With these 
conditions, the original networks are changed into feedforward networks without any feedback (reciprocial 
links or loops). Furthermore, the synchronization of the whole network is achieved in a hierarchical way. The 
conditions (i)-(iii) lead to the following path towards synchrony. There is a node in the directed spanning 
tree that has no input and acts as a master oscillator driving the network dynamics. If e {ai,a2), 
the oscillators that are just one hierarchical level below the master oscillator will synchronize. Then, the 
next lower level oscillators will also get synchronized and so on until the whole network reaches complete 
synchronization. Finally, we note that it could take a very long time to achieve complete synchronization 
when the number of effective layers is large. It would be interesting to see how the topology of the original 
networks is related to the depth of such effective directed trees and how it influences the transient time 
towards synchronization. Very recently, it was shown ^159l] that an age-based coupling similar to Eo. 11021 
but with Oij — e~"(*~^^^, will lead to such an effective directed tree with i? = 1 at large values of a [l59|. 

The findings of f96','l585 have important implications on the structural properties and dynamical processes 
in real networks. Although most analysis of complex networks have been developed for undirected networks, 
many real networks are directed. Recent studies about the local organization of directed networks found 



that reciprocity (percentage of bi-directional links) in real networks differ clearly from models [160|, Il6l| . 
and suprisingly many real directed networks have very few short loops as compared to random networks 
|l62| . These properties seem to be constrained by the degree correlations |l6lL Il62| , and therefore it will 
be interesting to study the impact of these properties on the synchronizability of directed networks in the 
future. 

Optimizing networks for synchronization has also been considered for networks of ph ase oscillators with 
natural frequencies uncorrelated with the initial random network configuration [g^, 163*1 . If the network is 
rewired with a bias towards connecting oscillators with similar average frequencies, synchronization is en- 
hanced and for intermediate coupling strength non trivial network properties, such as high number of cliques 
and large average distance, emerge 60^. When a measure that combines the local [39] and global synchro- 
nization is optimized by network rewiring, communities having similar natural frequencies are obtained at 
intermediate coupling, while strong coupling between dissimilar oscillators leads to highly synchronizale 
networks at large coupling strength |l63l |. 

Optimization approaches in networks thus are very useful to explore the structure-dynamics relation in 
model and real networked systems and also in various applications aiming at enhancing the performance of 



the networks, see |164| for a recent focus issue on this interesting topic. 
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Figure 29: The average values AX(/c) as a function of k in SF networks of Rossler oscillators outside the complete synchro- 
nization regime, (a) At various coupling strengths g = cr{k) below the threshold for complete synchronization and (b) the 
synchronized state is perturbed b y no ise of different intensities D. Qualitatively, the same behavior is observed when the 
oscillators are nonidentical. From III6II. 



4-3. Beyond the Master Stability Function formalism 

We have discussed how the impact of network topology on synchronizabiUty can be addressed using 
the MSF and graph theory. Away from the complete synchronization regime, the linear stability does not 
strictly apply. However, it is still possible to go one step forward to further understand some aspects of the 
dynamical synchronization patterns. 

In 11161] the authors studied effective synchronization patterns in unweighted SF networks of chaotic oscil- 
lators (with the coupling function H(x) = x) in several situations: away from the complete synchronization 
regime, when the coupling strength is smaller than the threshold for complete synchronization, when the 
oscillators have mismatches in parameters, and when there are noise perturbations. They considered a mean 
field approximation in which each oscillator is influenced by a global mean field X, with a coupling strength 
aki, namely, 

= F(x,) +f7fc»(X-x,), /c,>l. (106) 

The authors compared the synchronization of each oscillator to X by computing AXi = jx^ — X| and then 
obtained the average AX(fc) over all nodes with the same degree k. It was shown that out of the complete 
synchronized state 

AX{k) - k-'', (107) 

where the exponent 7 ~ 1, as seen in Fig. 1291 This result shows that in heterogenous networks, the hubs 
{ki > kth) will synchronize more closely with the mean field and they will form effective synchronization 
clusters (jx^ — Xj| < Ath)- However, there is not a unique threshold to define such effective clusters (see 
Fig. [^^. This path to synchronization, i.e., the formation of clusters by the hubs, was further described 
later in [3^, see also Sect. 13.1.41 

The formation of such SF or hierarchical clusters could be understood from a linear analysis using the 
MSF formalism. The linear variational equations of (jl06|) read 

i,^[DFp^)~ahI]^^, fc, »1. (108) 

They have the same form as except that Xi is replaced by ki and DH by the identity matrix I. The 
MSF for the couphng function H(x) = x is Amax(Q!) = Af — a, where Af is the largest Lyapunov exponent 
of the isolated oscillator F. Thus the largest Lyapunov exponent Aniax(fci) of the linearized Eq. ()108p is 
a function of ki and becomes negative for crki > Af . For large k values satisfying ak ^ Af , we have 
Amax(fc) « -o-fc. 
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Now suppose that the network is not completely synchronized, but slightly perturbed from the state of 
complete synchronization, when the coupling strength a is below the complete synchronization threshold, 
or when there is noise present in the system. For nodes with large degree fc, so that Amax(fc) ~ —crk is large 
enough in absolute value but negative, the dynamics of the averaged synchronization difference AX{k) over 
large time scales can be expressed as 



dt 



AX(fc) = A„,ax(fc)AX(fc)+C, 



(109) 



where c > is a constant denoting the level of perturbation with respect to the complete synchronization 
state, and depends either on the noise level D or on the coupling strength a. From this we get the asymptotic 
result AX{k) = c/|Amax(fc)|, yielding 

AX{k)^k-\ (110) 



which explains qualitatively the numerically observed scaling (solid lines in Fig. [29]). Interestingly, the same 
scaling dependence but for the time needed to get back into the fully synchronized state was obtained in 
[sit for a population of Kuramoto oscillators, see Sect. 13.1.41 

To round off this section, let us mention other works about s ynchronizability in complex networks that 
make use of linear criteria similar to the MSF 105 . 165 . Il66|. The analysis of the global stability of 



the synch ronized state, was first carried out for general graphs [l67 1 and then followed for specific comple: 



networks [iM llM llZO, llZll ■ The global stability requires additional constraints on the dynamical properties 
of the indivi dual osci llators. For example, consider the form of coupling as in Eq. (|54p . the requirement 
imposed by 168 . 169j | is that the following auxiliary system of the synchronization difference X^j = — Xj, 



D¥{pyij + (1 - f3)yii)dp - aDH 



be globally stable at the fixed point Xy 
synchronization of a network is 



for a > Uc- With this requirement, the condition for global 



a > (7 



max — 0/ , 



(111) 



where hi is the sum of the lengths of all chosen paths that pass through a given link I in the network. Note 
that bi is related to both the betweenness and the path length of the link. The result, which is for undirected 
networks, also holds for directed networks where the input and output degrees are equal for every node of 
the network. Finally, in 168 . 169l | it is also derived the condition needed for global synchronization in more 
general cases where each link in the network may have a different coupling strength that is allowed to vary 
in time. 



5. Applications 

The focus of the review up to now has been to revise the main contributions, from theoretical and 
computational points of view, to our understanding of synchronization processes in complex networks. In 
this section we will overview the applications to specific problems in such different scientific fields as biology 
and neuroscience, engineering and computer science, and economy and social sciences. There are nowadays 
several problems where the application of the ideas and techniques developed in relation to synchronization 
in complex networks are very clear and the results help to understand the interplay between topology and 
dynamics in very precise scenarios. There are other cases, also included here for completeness, for which 
most of the applications so far have been developed in simple patterns of interaction, but extension to 
complex topologies is necessary because it is its natural description. 

5.1. Biological systems and neuroscience 

In biology, complex networks are found at different spatio-temporal scales: from the molecular level up 
to the population level, passing through many intermediate scales of biological systems. In some of these 
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Figure 30: Dependency of the order parameter R on the the coupUng strength Q, which is linearly proportional to the cell 
density. The two curves correspond to different values of the ratio between mRNA and protein lifetimes variance A/3. From 




networks, dynamical interactions between units, which are crucial for our current understanding of living 
systems, can be analyzed in the framework of synchronization phenomena developed so far. Here we review 
some of these application scenarios where synchronization in networks has been shown to play an essential 
role. Thus, at the molecular level we can analyze the evolution of genetic networks and at the population 
level the dynamics of populations of species coupled through diffusion along spatial coordinates and through 
trophic interactions. Amongst these two extremes we find a clear application in the analysis of circadian 
rhythms. On a different context, neuroscience offers applications also at two different levels, one for the 
synchronization of individual spiking neurons and the other for the coupling between cortical areas in the 
brain. 



5.1.1. Genetic networks 

The finding that a few basic modules are the building blocks of large real regulatory networks has enabled 
the design and construction of small synthetic regulatory circuits to implement particular tasks. One of the 
most salient examples of a synthetic gene network is the "repressilator" , that has become one of the best 
studied model systems of this kind. The repressilator is a network of three genes, whose products (proteins) 
act as repressors of the transcription of each other in a cyclic way. This synthetic network was implemented 
experimentally in the bacterium E. coli, so that it p eriodically induces the synthesis of a green fluorescent 
protein as a readout of the repressilator state [172| . ft turns out that the temporal fluctuations in the 
concentration of each of the three components of the repressilator can be reproduced by a system of six 
ordinary differential equations, 



d[a 



dt ! + [%■]"' 

d[y. 



(112) 



dt 



- ~P{[y,]~[x,]), (lf3) 



where the couples (i, j) assume the values (1,3), (2, 1) and (3, 2). The variable [xi] is the mRNA concentration 
encoded by gene xi, and [yi] is the concentration of its translated protein j/^. The parameter a is the 
promoter rate, the parameter /3 is the ratio of the protein decay rate to the mRNA decay rate, and time has 
been rescaled in units of the mRNA lifetime. This system has a unique steady state which can be stable or 
unstable depending on the parameter values and constitutes an illustrative example of the experience gained 
by identifying network modules and modeling its dynamical behavior in real networks. 

Not surprisingly, the repressilator has called much attention from experts on biological synchronization, 
since it offers good prospectives for further insights into the nature of diverse biological rhythms -whose 
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mechanisms remain to b e understood- which are generated by thousands of ceUular oscillators that operate 
synchronously. In |l73| it has been recently proposed a simple modular addition of two proteins to the 
repressilator original design that could be used to describe the metabolic oscillations observed in a well 
mixed suspension of yeast cells. In the new setting, one of the new proteins can diffuse through the cell 
membrane, thus providing a coupling mechanism between cells containing repressilator networks. This 
inter-cell communication couples the dynamics of the different cell oscillators (with different repressilator 
periods) and thus allows us to study the transition to synchronization of coupled phase oscillators in a 
biological system. In particular, in the limit of infinite cell dilution, the system is made up of a population 
of uncoupled limit-cycle oscillators. This is not anymore the case when the cell density increases, as now 
extracellular diffusion provides a mechanism of intcrcell coupling. As a result, the system shows partial 
frequency locking of the cells. Finally, if the cell density is large enough, complete locking and synchronized 
oscillations are observed. The authors proposed an order parameter to measure the degree of synchronization 
of oscillatory behavior. The dependence of this order parameter R defined as the ratio of the standard 
deviation of the time series of the average signal to the standard deviation of each individual signal [xi] 
averaged over all signals i, as a function of the coupling strength, is shown in Fig. [30] for different values of 
the ratio between the mRNA and protein lifetimes width distribution {6/3). Note that the phase diagram 
of the coupled repressilators can be explained by the very same mechanism involved in the transition to 
synchronization for systems of coupled oscillators (e.g., Kuramoto oscillators) studied in Sec. 13.11 What 
is relevant here is that the transition from an unsynchronized to a synchronized regime is caused by an 
increase in cell density and therefore the experimental observation of a synchronizing transition in biological 
phase oscillators might be achieved. In fact, it has b een recently designed a simple electronic circuit analogy 



of a population of globally coupled repressilators [174|. They show that coupling is more efficient than 



externally forcing for the achieveme nt of synchronization. In contrast to the existence of a unified rhythm 
that gives rise to synchronization, in [l75| it is analyzed the mechanisms of intercell communication that can 
be responsible for multirhythmicity in coupled genetic units. We foresee that works on this line of research 
will incorporate more genetical interactions in the near future, being the complex network substrate and the 
synchronization dynamics key aspects of the whole problem. 

5.1.2. Circadian rhythms 

A circadian rhythm is a roughly 24-hour cycle in the physiological processes of living systems; usually 
endogenous, or when it is exogenous it is mainly driven by daylight. Understanding circadian rhythms is 
crucial for some physiological and psychological disorders. A nice description of experim ents carried o ut in 
human beings, in which their circadian rhythms are altered, can be found in the books bv lStrogat j 17d | and 



Pikovskv et al.l [l[. Circadian rhythms are known to be dependent on the network of interactions between 



different subsystems. For example, daylight sensed by eyes and processed by the brain develo ps a chain of 
interactions that affects even the behavior of certain groups of cells. On a different scenario, jl77l| reports 
how non-oscillatory cardiac conducting tissues, when driven rhythmically by repetitive stimuli from their 
surroundings, produce temporal patterns including phase locking, period-doubling bifurcation and irregular 
activity. 

Synchronization phenomena in complex networks of coupled circadian oscillators has been recently inves- 
tigated experimentally [l78l | on plant leaves. The vein system is in this case the complex network substrate 
of the synchronization process. Plant cells are coupled via the diffusion of materials along two types of con- 
nections: one type that directly connects nearest-neighboring cells and the other type that spreads over the 
whole plant to transport material among all tissues quickly. Analyzing the phase of circadian oscillations, 
the phase-wave propagations and the phase delay caused by the vein network, the authors describe how 
global synchronization of circadian oscillators in the leaf can be attained. As we have seen throughout this 
review, the role of the topology of interactions is again fundamental in the development of synchronization. 
This work is representative of the new type of applications we can find in the very recent literature about 
synchronization in complex networks. This particular case of circadian rhythms in plants might be extended 
to other living systems, including humans. 
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5.1.3. Ecology 

It is a well known observation in nature that fluctuations in animal and plant populations display complex 
dynamics. Mainly irregular, but some of them ca n sho w a remarkably cyclical behavior and take place over 
vast geographical areas in a synchronized manner |l79l | . One of the best documented cases of such situation 
are the population fluctuations in the Canadian lynx, obtained from the records of the fur trade between 
1821 and 1939 in Canada. Fluctuations in lynx populations show a 10-year periodic behavior from three 
different regions in Canada On the other hand, there are some evidences that the existence of 



conservation corridors favoring t he dis persal of species and enhancing the synchronization over time increases 
the danger of global extinctions 

One of the first explanations for such types of behavior was that of synchronous environmental fo rcing , 



this is the so-called Moran effect 184.] • There are, however, other explanations for this phenomenon [185|, 



but in any case the problem highlights the impor tance of integrating explicitly spatial and trophic couplings 



into current metacommunity theories [186l |187[. Some efforts along these lines have already been made 
by c onsidering very simple trophic interaction in spatially extended systems. For example, the authors in 
[181 ' analyze a three-level system (vegetation, herbivores, and predators), where diffusive migration between 
neighboring patches is taken into account. They find that small amounts of m igration are required to induce 
broad-scale synchronization. Another interesting study is performed in [l88j . where, again in an extremely 
simple model, it is found that changing the patterns of interaction between consumers and resources can 
lead to either in-phase synchrony or antiphase synchrony. 

Nowadays we know, however, about the inherent complexity of food- webs jl89l |. Food webs have been 
studied as paradigmatic examples of complex networks, because they show many of their non-trivial topo- 
logical features. Furthermore, the existence of conservation corridors affecting the migration between regions 
adds another ingredient to the structure of the spatial pattern. It is precisely this complexity in the trophic 
interactions coupled to the spatial dependence that must to be considered in detail in the future to get a 
deeper understanding of ecological evolution. 

5.1.4. Neuronal networks 

Synchronization has been shown to be of special relevance in neural systems. The brain is composed of 
billions of neurons coupled in a hierarchy of complex network connectivity. The first issue concerns neural 
networks at the cellular level. In the last years, significant progress has been made in the studies about 
the detailed interconnections of different types of neurons at the level of cellular circuits T90, 191, 192]. At 
this level, the neuronal networks of mammalian cortex also possess complex structure, sharing SW and SF 
features. Here are two basic neuron types: excitatory principal cells and inhibitory interneurons. In contrast 
to the more homogeneous principal cell population, interneurons are very diverse in terms of morphology and 
function. There is an apparently inverse relationship between the number of neurons in various interneuron 
classes and the spatial extent of their axon trees-most of the neurons have only local connections, while 
a small number of neurons have long-range axons [193] . These properties of neuronal networks reflects a 
compromise between computational needs and wiring economy jl94l . Il95| . 

On the one hand, the establishment and maintenance of neuronal connections require a significant 
metabolic cost that should be reduced, and consequently the wiring length should be globally minimized. 
Indeed, the wiring economy is apparent in the distributions of projection length in neural systems, which 
show that most neuronal projections are short 196l . 197 1. However, there also exists a significant number of 



snow tnat most neuronal proiectioni 
long-distance projections [l98l . Il08j . 



Large-scale synchronization of osc illat ory neural activity has been believed to play a crucial role in the 
information and cognitive processing jl99l |. At the level of cellular circuits, oscillatory timing can transform 
unconnected p rinci pal cell groups into temporal coalitions, providing maximal flexibility and economic use 
of th eir spikes [20Cll | . Brains have developed mechanisms for keeping time by inhibitory interneuron networks 



20l| . The wiring will be the most economic if the connections were all local. However, in this case physically 
distant neurons are not connected, and synaptic path length and synaptic delays become exceedingly long 
for synchronization in large networks. From previous analysis of synchronization of random networks, we 
know that synchronizability (stability of the synchronized state) i s op timal in fully random networks with 



a uniform connectivity per node, independent of the network size lllj . The same happens if interneuronal 
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oscillators are coupled [193|. However, fully random connections irrespective of physical distance are not 
econo mic if wiring cost is taken into account. 

In [l93|. it was shown with a model of interneuronal networks containing local neurons (Gaussian dis- 
tribution of projection length) and a fraction of long-range neurons (power law distribution of projection 
length) , that the ratio of synchrony to wiring length is optimized in the SW regime with a small fraction of 
long-range neurons. Thus, most wiring is local and neurons with long-range connectivity and large global 
impact are rare, as consistent with observations. It was argued that the complex wiring of diverse interneu- 
ron classes could represent an economic solu tion for supporting global synchrony and oscillations at multiple 
time scales with minimum axon length [l93l | . While such mathematical consideration can predict the scaling 
relationship among the interneuron classes in brain structures of varying sizes, understanding the role of 
complex neuronal connectivity, most likely mediated by synchronization, is still one of the main challenges 
in neuroscience. The theory reviewed in this article will surely contribute to their understanding when more 
systematic information of neuronal connectivity becomes available in the future. 



5.1.5. Cortical networks 

The application of graph theoretical approaches, and the characterization of the functional activity by 
neural sy nchronization, have significantly contributed to the understanding of complex networks in the brain 



On a larger neurophysiological scale, the activity observed experimentally by electroencephalographs or 
functional magnetic resonance imaging, is characterized by oscillations occurring over a broad spectrum 
and by synchronization phenomena over a wide range of spatial and temporal scales. Reliable databases are 
available now for large-scale sys tems leve l connectivity formed by long-range projections among cortical areas 
in the brains of several animals 142l 204 1. Large-scale b rain networks a re found to be densely connected, with 



very complex and heterogeneous connectivity patterns [2051 . 12061 . 1207f . In parallel, the investigation of brain 
activity has also put significant emphasis on large-scale funct iona l interactions, c haracterized by coherence 
and synchronization between the activity of cortical regions 208 . 209L 210L 211 1. Both the structural and 
functional connectivity of the brain display SW and SF features. The relationship between structural and 
functional connectivity remains an important open problem in neuroscience. 

Recent simulations of sync hronization dy namics of brain networks have shed light on this challenging 
problem. In a series of papers [212l I213L l214l | the dynamics o f a re alistic cortico-cortical projection network 
of the cat has been modeled at the level of funct ional areas [l42l |. At this level, the network (see Fig. [31]) 
displays a hierarchical cluster organization |207| . There are four prominent clusters that agree broadly 
with the four functional subsystems: visual (V), auditory (A), somatomotor (SM), and frontolimbic (FL). 
They simulated the network dynamics by a 2-level model: each node (cortical area) is represented by a SW 
subnetwork of neurons (network of networks). It was shown that the model possesses two distinguished 
regimes, weak and strong synchronization. In the weak synchronization regime, the model displays biologi- 
cally plausible dynamical clusters. The functional connectivity, obtained by passing the correlation matrix 
through various thresholds, exhibits various levels of organization. The clusters with the highest levels of 
synchronization are from respective functional subsystems (Fig. [35] (a,b)) and are related to specialized 
functions of the subsystems. The specialized clusters are integrated into larger clusters through brain areas 
having many inter-community connections (Fig. [32| (c,d)). As a whole, the functional connectivity reveals 
the hierarchical organization of the structural connectivity [212| . The dyn amics forms four major clusters 
(Fig. [33]), in excellent agreement with the four functional subsystems 213[. Furthermore, brain areas that 
bridge different dy nami cal clusters are found to be the areas involved in multisensory associations. In a 
comparative study |214l |. it was shown that repr esen ting the brain areas with a periodic, low-dimensional 
neuronal mass oscillator describing alpha waves [215l | cannot resolve these four clusters. The detailed net- 
work topology becomes rather irrelevant to the dynamical patterns which is not very much changed when 
the network is randomized. This is t he sa me for the strong synchronization regime in the 2-level model 
which resembles epileptic- like activity |213| .). This can be understood recalling previous analysis based on 



random networks where is shown that synchronization is mainly determined by the node intensity 111], |117 1 



(see Section IV 6). Furthermore, the transition from the weak to the strong synchronization regime shares 



a similar picture with the Kuramoto model in complex networks [39|,]40| (see Section III 4). 
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Figure 31: (Color online) Connection matrix AI^ of the cortical network of the cat brain. The different symbols represent 
different connection weights: 1 (o sparse), 2 (• intermediate) and 3 (* dense). The organization of th e system into four 
topological communities (functional sub-systems, V, A, SM, FL) is indicated by the solid lines. From j212| '). 
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Figure 32: (Color online) The functional networks (o) at various thresholds: Rth = 0.07 (a), Rt^ = 0.065 (b), Rt^ = 0.055 
(c) and Rth = 0.019 (d). The small dots indicate the anatomical connections. From |212|| . 
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Figure 33: (Color online) Four m ajor dynamical clusters (o) in the weak synchronization regime, compared to the underlying 
anatomical connections (■). From |213j . 



Shortly after these works, a very similar structure-function relationship was observed [216|. Each area of 
the macaque neocortex was represented by a neural mass model in the regime of spontaneous activity with 
complicated temporal patterns. It was shown that the functional connectivity, measured over a very long 
time, is closely shaped by the underlying structural connectivity as described in 212 . 213L 214l |. On short 
time-scales, the functional connectivity changes, forming two anticorrelated clusters, similar to functional 
networks obtained from brain imaging data (217.) . 

These findings support the idea that the brain is an active network, and it can generate activity by 
itself in the absence of external signals. Classical theories in cognitive neuroscience viewed the brain as a 
passive, stimulus- drive n device and the spontaneous on-going activity of the brain had been regarded as 
background noise [218| . It is still customary in data analysis to take the average signals over many trials of 
electroencephalographs as the event-related activation and associate them with cognitive processes. In the 
view of active dynamical brain networks, it has been sh own that the spontaneous on-going activity imposes 
significant impact on the selective responses to stimuli 218L 199l |. The intricate re l ation ship between large- 

[ili . [ili p will contribute to 



212, m 



scale structural and functional networks revealed in these works 
this reorientation of concepts in neuroscience. On the other hand, excessive activation and synchronization 
of neural networks have been found to associate with dysf unctions and disorders of the brain, such as 
the epileptic seizure [219| and the Parkinsonian disease [220l |. Understanding synchronization in neuronal 
networks of various level, especially studying the role played by the complex network topology, is crucial 
to elucidate how normal brains can maintain desirable levels of synch roniz ation. It will also contribute 
significantly to biomedical data analysis of pathological brain activities 22l| . for example, the challenging 
task of detecting precursors that can make prediction much before the clear onset of seizures , and design 
suitable methods for treatments of neural diseases 12221. 



5.2. Computer science and engineering 

Complex networks and synchronization dynamics are relevant in many computer science and engineering 
problems. For example, in computer science, synchronization is desirable for an efficient performance of 
distributed systems. Sometimes, the goal of the distributed system is to achieve a global common state 
(consensus). Nowadays these systems are becoming larger and larger and their topologies more and more 
complex. On the other hand, some engineering problems also face the need of maintaining coordination at 
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Figure 34: (Color online) Virtual time horizon snapshots in the steady state for 10000 sites in one dimension, (a) For a regular 
network. The lateral correlation length g and width to are shown for illustration in the graph, (b) For the SW n etwo rk with 
p = 0.1, the heights are effectively decorrelated and both the correlation length and the width are reduced. From |224l . 



the level of large scale complex networks, for example in problems of distribution of information, energy or 
materials. 



5.2.1. Parallel/ Distributed computation 

The simulation of large systems are, nowadays, mainly implemented as parallel distributed simulations 
where parts of the system are allocated and simulated on different processors. A large class of interacting 
systems including financial markets, epidemic spreading, traffic, and dynamics of physical systems in general, 
can be described by a set of local state variables allowing a finite number of possible values. As the system 
evolves in time, the values of the local state variables change at discrete instants, either synchronously or 
asynchronously, depending on the dynamics of the system. The instantaneous changes in the local configura- 
tion are called discrete events, forming what has been coined as a parallel discrete-event simulation (PDES) 
[223| . The main difficulty of PDES is the absence of a global pacemaker when dealing with asynchronous 
updates. This imposes serious problems because causality and reproducibility of experimental results are 
desired. In a conservative scheme, processes modeling physical entities are connected via channels that 
represent physical links in the target system. Since PDES events are not synchronized via a global clock, 
they must synchronize by communication between nodes. 

The essentials of a PDES are: a set of local simulated times of the processors and a synchronization 
schedule. As the grid-computing networks with millions of processors emerge, fundamental questions of the 
scalability of the underlying algorithms must be addressed. The properties of a PDES to be scalable are: 
(i) the virtual time horizon {Ti(i)}f'', a set of time simulated instants in Np processors after t parallel steps, 
should progress on average with a non-zero rate, and (ii) the width of the time horizon should be bounded 
when N„ ^ oo. 



In j225l | , it was studied the asymptotic scaling properties of a conservative synchronization algorithm in 
massive PDESs where the discrete events are Poisson arrivals. They found an interesting analogy between 
the evolution of the simulated time horizon and the growth of a nonequilibrium sur faceS They showed that 
the steady-state behavior of the macroscopic land scap e of the simulated time horizon, in one dimension, is 
governed by the Edwards- Wilkinson Hamiltonian 228l | . 



*Note th at th e first analogy between synchronization processes and the theory of surface growth appeared in |226|| , posteriorly 
revisited in [227l | (see [l[ for a comprehensive review). 
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The analogy becomes clear by interpreting the virtual times as the height of a surface, and defining 
the width of the simulated times surface as the mean root square of the virtual times with respect to the 
mean r. This width provides a measure for de-synchronization 




(114) 



The problem that now is faced is how the width of the synchronization landscape scales with the number 
of processors. If the scaling diverges, it means that the synchronization is hardly achievable. In one- and 
two-dimensional regular lattices, the width of the synchronization landscape diverges with the number of 
processors as w ~ nI^'^ where d is the dimension. This effect can be traced b ack to the lateral correlation 
length ^ of the surface that also diverges with the numb er of processors ^ N [229| . An interesting solution 
to this problem has been proposed in [230l | and [224| . It consists of adding a few random links to the 
regular lattices resulting in a SW structure. This structure has the effect of de-correlating the lateral length, 
suppressing large fluctuations in the synchronization surface (roughness), and producing finite average values 
of w in the large system-size limit, see Fig. 1341 Moreover, the extreme height diverges only logarithmically in 
this limit. This latter property ensures synchronization in a practical sense in a SW topology of processors. 



5.2.2. Data mining 

The term data mining refers to the process of automatical ly se arching large volumes of data for patterns 



that provide some useful information for classification. In 23lj it has been proposed a new method of 
data mining based on spontaneous data clustering using a network of locally coupled limit-cycle phase 
oscillators. The method is closely related to the determination of community structure via synchronization 
processes devised by several authors [H^ |6lj . The idea is to encode multivariate data vectors (that are the 
elements of the database) into vectors of natural frequencies for an oscillators' dynamical model, akin to 
the KM, expecting that the dynamics of the system will group similar data in clusters of synchronization. 
More precisely, given n multivariate data points with m degrees of freedom, xl — {xi{l), Xi{2), Xi{m)), 
i ~ 1, . . . , n, they write the dynamical model: 

n 

^,(0 =x,(0 + ^5I^(^'J>i^(^j(0-^.(0) (115) 

where 9i{l) is the l-th component of the phase vector 9i = {9i{l),9i{2), ...,9i{m)), H{dij) is a function 
that determines the neighborhood of 9i based on the distance dij — \xi — Xj\. The determination of the 
neighborhood provides the network of interactions between oscillators. The proposal by the authors is a 
neighborhood centered at Xi defined by the hyper-sphere of radius dg, being do — a\xi\ and a a tuning 
parameter. The function H(dij) = 1 if dij < dn and otherwise. The application of the method in a 
database of aging status in frail elderly reported in |23l| shows a good performance of the method, and gives 
a nice expectative of exploitation of the concepts of synchronization in the area of data mining. 



5.2.3. Consensus problems 

Consensus problems, understood as the ability of an ensemble of dynamic agents to reach a unique 
and common value in an asymptotically stable stationary state, have a long history in the field of computer 
science, particularly in automata theory and distributed computation. In many applications, like for instance 
cooperative control on unmanned air vehicles, formatio n con trol or distributed sensor networks, groups of 
agents need to agree upon certain quantities of interest |232| . As a result, it is important to address these 
problems of agreement within the assumption that agents form a complex pattern of interactions. These 
interactions can be directed or undirected, fixed or mobile, constant or weighted, keeping then many of the 
ingredients we have been discussing in this review. Another interesting fact in this sort of problems is the 
existence of time delays in the communication process. 
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In |232l | the authors define consensus problems on general graphs. Let us consider a dynamic graph in 
which the connectivity pattern of the nodes can change in time. At each node, a dynamical agent evolves 
in time according to the dynamics 

Xi^f{xi,Ui), (116) 

where f(xi,Ui) is a function that depends on the state of the unit Xi, and on Ui that describes the influence 
from the neighbors. The x-consensus problem in a dynamical graph is a distributed way to reach an 
asymptotically stable equilibrium x* satisfying x* = x(x(0)), Vz, where x(2;(0)) is a prescribed function of 
the initial values (e.g., the average or the minimum values). 

They present two protocols that solve consensus problems in a network of agents: 

1. fixed or switching topology and zero communication time-delay: 

N 

Xi=^ aij{t){xj{t) - Xi{t)), (117) 

2. fixed topology and non-zero communication time-delay Tij > 

TV 

= ^ aij{xj{t - Ty) - Xi(t- Tij)). (118) 

ij = l 

We note that the analysis of the asymptotic behavior of such linear system is similar to the stability analysis 
performed in the framework of the MSF (Sect. [4]). 

The authors find very interesting results in terms of network properties. For instance, a network with 
fixed topology that is a strongly connected digraph (a subgraph connected via a path that follows the 
direction of the edges of the graph), solves the average consensus problem if and only if all the nodes of the 
network have the same indegree as the outdegree, i.e. fcj" = fc™*, Vz, as the balanced networks discussed in 
[l69l | (see Sect. 14. 3p . Furthermore, the performance of the network measured in terms of the speed in which 
the global asymptotic equilibrium state is reached, is proportional to X2{Q) where Q is the mirror graph 
induced by Q, which is defined as an undirected graph with symmetric adjacency matrix dij = {flij -\- aji)l2. 

For a switching topology, they find that if the dynamics of the network is such that any graph along the 
time evolution is strongly connected and balanced then the switching system asymptotically converges to an 
average consensus. Concerning time communication-delays, the important result is that if all links have the 
same time-delay r > 0, and the network is fixed, undirected and connected, the system solves the average 
consensus if r e (0,7r/2Ajv). In this case, in a similar way as discussed in previous applications, there are 
two tradeoff issues that can be related to problems of network design; one concerns the robustness of the 
protocol with respect to time-delays, and the other to communication cost. 

When applying this framework to a certain class of networks, it is found | , 233l| that the speed of conver- 
gence, as the inverse of A2 (as was also found in synchronization problems [52]), can be increased by orders 
of magnitude by simply rewiring a regular lattice, while this change has a negligible effect on Aat, which 
measures the robustness to delays of the system. This can be understood by the eigenvalues of the SW 
network in Eq. (|70p as compare d to the regular networks in Eq. (j67p . Some other variations can be found 



in the recent literature; e.g., in 234j several network models with physical neighborhood connectivity are 
analyzed. Depending on the precise rules they discuss on the performance and the robustness of the system. 

Due to its importance as an application in computer science, consensus problems are interesting by 
themselves. But understanding its linear dynamics can be also of great importance in the behavior of 
complex populations of units that evolve according to more complex non-linear dynamics, as it happens in 
many synchronization problems. 

5.2.4- Wireless communication networks 

Another emerging line of research can be found in the field of synchronizing wireless sensor communication 
networks. Wireless ad-hoc networks are telecommunication infrastructures formed by devices equipped with 
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a short-range wireless technology, such as WiFi or Bluetooth. Unlike wired networks, these network s can 



be created on the fly to perform a task, such as information routing, environmental sensing, etc. 235 |. 
Furthermore, the topology of these networks can be changed dynamically to achieve a desired functionality. 
From the perspective of fundamental research, these systems provide a clear-cut example of highly dynamic, 
self-organizing complex systems. 

One of the main technological problems in wireless networks is that of synchronize time activity in a 
decentralized way. Wireless time synchronization is used for many different purposes including location, 
proximity, energy efficiency, and mobility for example. We will revise in this section two approaches to 
the problem that have been developed within the scope of the synchronization scenario reviewed so far. 
Other approaches not clearly connected to synch roni zation in complex networks to solve this important 
problem can be found in the specific literature [236| . The first approach concerns to the routing and 
information flow algorithms which require synchronization of the clocks of the nodes of the wireless network 
to establish a global coordinated time. The topology of these sensor networks is accurately represented by 
random geometric graphs which are constructed by dropping n points randomly uniformly into the unit 
square (or more generally according to some arbitrary specified density function on d-dimensional Euclidean 
space ) and adding edges to connect any two points distant at most r from each other. In a very recent 



work [237[, synchronization of Kuramoto oscillators in these networks is studied. They consider a wireless 



system in which the connections vary at a time scale much shorter than the time scale associated to the 
synchronization dynamics, and hence the network is static. Nodes correspond to devices that have a finite 
transmission range, and are linked to those nodes that are located within the range. This procedure gives 
rise to a two-dimensional random geometric graph, which is characterized by a high clustering coefficient 
and a very large average shortest path length, as compared to ER graphs with the same number of nodes 
and links. The remarkable result is that this type of network is very hard to synchronize, both in terms of 
the stability of the synchronized state and in terms of the time needed to reach the completely synchronized 
state. Although they are very homogeneous, the smallest non-zero eigenvalue of the Laplacian matrix is 
very low, providing a clear limitation for synchronizability of Type II. This result points out the limitations 
concerning synchronizability that raise from this topology, interestingly, just by rewiring a small fraction of 
the links at random synchronizability is clearly improved, the distances are shortened and at the same time 
the clustering is decreased, which show up as an increasing of the eigenvalue A2 almost without affecting the 
largest eigenvalue Aat. An extended study about the convergence properties of a similar s ystem where nodes 
are represented as discrete-time oscillators, is studied through algebraic graph theory in [238j. The authors 
main finding is that the distributed synchronization problem converges to a unique cluster of synchronized 
nodes, if and only if the associated weighted directed graph G is strongly connected, i.e. if there is a path 
from each vertex in the graph to every other vertex. 

The second approach concerns the shared resources available in these systems. Communication channels 
have a finite bandwidth so that the access times of different users should be desynchronized when their 
number is large and non necessarily constant, this is basically the idea behind any Time Division Multiple 
Access (TDMA) algorithm to be used over a multi-hop wireless network. In [239|, the authors proposed 
a biologically inspired algorithm for desynchronization in a single-hop network that is in the scope of the 
review. They consider a set of N nodes (integrate-and-fire oscillators) that generate events with a common 
period. The nodes rearrange their phases, just considering the times in which neighboring (in time) units 
fire, in such a way that the events become spaced at intervals T/N . The final state then corresponds to what 
is usually called a round-robin schedule. In this w ay, th e use of the bandwidth without collisions between 
messages is optimized. Inspired by this result, in |240 | the authors considered the units to be Kuramoto 
oscillators with a common frequency. Introducing some dephase angle in the sinus function and coupling 
pairs of units along a closed chain, the authors find new stable configurations different from the completely 
synchronized state. Some of these configurations correspond also to the round-robin schedule, which turns 
out to be very robust under the addition or deletion of nodes from the system. This finding obtained in the 
scope of a precise application, is general, and accounts for global synchronization of discrete-time dynamical 
directed networks. 

To end this section we also to point out a different problem, also in the area of wireless mobile sensors, 
that has been nicely solved by a description in terms of pulse coupled oscillators in complex networks. The 
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problem is that of the decentralized detection of abrupt changes, that in wireless networks can represent 
failures in communication, attacks to certain sensors or, generally speaking, any change in the sensed activity 
of purpose. The authors of [241*1 proposed a very simple transmitter with no routing, no multiple access, 
only a "pulse position modulation" mechanism. In particular, they assume that the nodes can transmit 
only through the emission of pulses with constant amplitude. The information of the sensor data and 
the interaction among them can only be encoded in the timing of the pulse emission. This simplistic but 
effective configuration leads the problem of decentralized detection of abrupt changes to that of observing 
the synchronization of pulse coupled oscillators in a network, after a local perturbation of them. The work 
did not included the specific topology of interactions as a fundamental parameter of the problem, and then 
do not propose an optimum network scenario for the propagation of the signals, however it clearly points 
out another technological problem in electrical and computer engineering, that can be solved in the scope 
of synchronization in complex networks. 

5.2.5. Decentralized logistics 

Logistics is the management of the flow of goods, information and other resources. Sometimes, this 
management is limited by the capability of maintaining global knowledge and/or global communication. In 
this case, the necessity of decentralized coordination mechanisms are mandatory. In many material flow 
systems coordination of tasks in a parallel way is essen tial f or an optimal functioning but difficult to ac hieve . 
Typical examples of this are cross-roads in road traffic [242l | and supply chains in production processes |243l | . 




Recent work on supply networks has shown how to treat them as physical transport problems governed by 
balance equations and equations for the adaptation of production speeds. Although the nonlinear behavior 
is different, the linearized set of coupled differential equations is formally related to those of mechanical or 
electrical oscillator networks [244] . Whereas traditional optimization techniques can be used to setup single 
nodes, the inherent topological complexity makes maintenance of coordination at network-wide level to be 
practically unsolvable by these methods. Furthermore, robustness and ffexibility, due to continuous changes 
in demand and failures, are also required for an efficient transportation. 

There is an analogy between material transportation in networks and the flow of chemical substances 
and nutrients in biological organisms, where synchronization dynamics plays an important role. In 24^ it 
has been proposed a decentralized control model for material flow networks with transportation delays and 
setup-times, based on phase-synchronization of oscillatory services at the network nodes. 

A material transport network is a directed and weighted graph where the flow of material at nodes is 
conserved. Subsets of links are active at different times, and this makes that the activity of the node is 
periodic and one can map a continuous phase variable 0{t) to a discrete service state M : 9 it) — > s{t). While 
the phase angle 6 of the oscillator modeling the intersection varies from to 2tt at a rate w, all states s 
are sequentially activated. To achieve a coordination of the switching states on a network-wide level, they 
propose a suitable synchronization mechanism. 

The authors apply this formalism to the control of traffic lights at intersections of road networks. A 
single traffic light intersection is modeled by an oscillator where the continuous phase maps to a set of states 
corresponding to green lights (see Fig. I35p . The maximum frequency of the oscillator dynamics is calculated 
in terms of the load at the different lanes and the setup-time. Global coordination of the network is achieved 
by synchronizing the local phases and frequencies, requiring to reach a phase-locked state where the phase 
difference between neighboring sites is fixed. They suggest a coupling on two different time scales: 

• Adaptation of the phase 6* a la Kuramoto: 



where oji is the intrinsic frequency. As long as u)i < ojf^'^^, 6i tries to adjust to the neighboring phases. 
The constant Tg corresponds to the typical time scale for this adaptation. 





(119) 
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Figure 35: (Color online) A single intersection adjusts the mapping of the phase-angle 6 to the switching states s locally. Within 
a complete cycle, each state s is sequentially activated for a period AOs, during which the corresponding non-conflicting traffic 
lights are set to green. While switching from one state to another , all traffic lights are set to red for a period of A9^^*^P. The 
phase-angle, at which a new cycle starts, is denoted by 9o. From [245j . 



• A second decentralized coupling can be used to increase the intrinsic frequencies to approach the 
possible maximum within a slow time scale: 

uj^ = ^ (mm{ujj{t)} + Auj~ uj,{t)] . (120) 

Here the constant parameter Alo provides a linear drift towards higher frequencies. 

Under these assumptions two dynamical regimes are possible (see Fig. I36|) . Starting with a random 
initial condition (left), the system quickly settles into a state with growing common frequency and vanishing 
phase-differences. As soon as the maximum common frequency is found, the system enters the other state 
with a locked common frequency and phase-differences exponentially converging towards constant values 
(right). 

The extension of the analysis performed on su ch a simple setting to more realistic complex transportation 
networks is very promising but challenging |24d |. 

5.2.6. Power-grids 

Power grids are physical networks of elec trical power distribution lines of generators and consumers. In 
the pioneering paper by .Watts and Strogat3 ^ it was already reported that the power-grid constitutes one 
of the examples of a self-organized topology that has grown without a clear central controller. This topology 
is indeed very sensitive to attacks and failures. From its topological point of view there are several analyses 
on power-grids in different areas of the world 2471 . |248| | and some models have been proposed to deal with 



the cascading process of failures EH^ ^EM- 



The principles of electricity generation and distribution are well known. Synchronization of the system 
is understood as every station and every piece of equipment running on the same clock, which is crucial for 
its p roper operation. Cascading failures related to de-synchronization can lead to massive power blackouts 
[25^ . 

Then power production, dissipation, transmission, and consumptio n rep resents a dynamical problem and 
the power grid can be seen as an example of a system of oscillators [l76l |. Along this line, in [253l | it has 
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Figure 36: A single intersection of roads is modelled by a phase angle 9. This continuous phase maps to a discrete set of states 
{s\. Simulation results for a regular lattice road network, where the intersections are defined as oscillators with (a) a frequency 
u)i and (b) a phased;, shows a phase transition towards a synchronized (frequency-locked) state. From 12451 . 



beeb proposed a model where each element (generators and machines) is described in terms of a phase that 
grows linearly in time with a frequency that is close to the standard frequency of the electric system (50 or 
60 Hz). 

Consider the power produced at a generator. It can then be dissipated, accumulated, or transmitted 
along the electric line. The first two terms (dissipation and accumulation) depend on the frequency of the 
generator whereas the last one (transmission) is proportional to the sinus of the phase difference between 
the generator and the machine at the other extreme of the line. Then, a simple energy balance equation 
relates the evolution of the phase (first and second time derivatives) with sinus of phase differences. 

Applying this simplified approach to a networked system of generators and machines, they arrive to a 
set of Kuramoto-like differential equations 

e, = -ae, + LU, + K J2 sin(ej - 0,), (121) 

where LUi is related to the power generated at the element and to the dissipated power, and K, representing 
the stength of the coupling, is related to the maximum transmitted power. 

Within this framework, they analyze, as an application, under which conditions the system is able to 
restore to a stable operation after a perturbation in simple networks of machines and generators. To 
the best of our knowledge this is a first approximation to the real applicability of the knowledge about 
synchronization in complex networks to power grids, although the hypothesis along the work are still very 
relevant here. 

5.3. Social sciences and economy 

In the last decades, social sciences and economy have become one of the favorite applications for physi- 
cists. In particular, t ools and models from statistical physics can be implemented on what some people has 
called social atoms, j254l | i.e. unanimated particles are replaced by agents that take decisions, trade stocks 
or play games. Simple rules lead to interesting collective behaviors and synchronization is one of them, 
because some of the activities that individual agents do can become correlated in time due to its interaction 
pattern, which, in turn, is clearly another example of the complex topologies considered along the review. 

In social systems, however, it is not an easy task to identify the relationship between agents (being 
humans or collectives in social interactions, stock prices in finances, or countries in the World Trade Web). 
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We alert the reader, that in some cases the quasi-periodic correlated activity is interpreted as synchronization 
in this scenario. For example, dynamical similarity along economic cycles is understood as synchronization 
in economic terms. Keeping this in mind, some of the applications presented in this section are weak 
formulations of synchronization in complex networks, however we think they are interesting because they 
open the door to stronger formulations in this context. 



5.3.1. Opinion formation 

One of such problems is the study of opinion formation in society. The underlying idea is that individuals 
(or agents) have opinions that change under the influence of other individuals giving rise to a sort of collective 
behavior, grouping together a macroscopic part of the whole population with similar social features [l5l |. 
Therefore, the main goal is to figure out whether and when a complete or partial consensus can emerge out 
of initially different opinions, no matter how long it takes for the consensus to be reached. 

In general, the formation of a collective opinion about a certain matter is not equivalent to a transition 
to som e kin d of synchrony, but rather to a transition to an absorbing state. However, a recently proposed 
model |255l | makes explicit use of a modified KM (see Sect. Ill) and thus in this case the formation of 
groups of opinions can be thought of a synchronization process. Agreement models deal with N individuals 
characterized by an opinion Xi (either an integer or real number) and a network of contacts that drives 



the dynamics of opinion formation through deterministic rules (15| . In 255l | it has been considered the 
case in which opinions are neither bounded nor periodic, but that two initially different opinions can also 
diverge when time goes on. Moreover, they also take into account that two quite different individuals tend 
to interact less by assuming that the coupling between these two individuals is a decreasing function of their 



opinion differences. Finally, in 255l | the main source of heterogeneity is not given by the initial positions, 
but by different rates of changing individuals' opinions. 

Taking all the previous statements into account, it has been proposed the following governing equations 
for the dynamics of the rate of change of opinions ^255] 

N 

iz = c^. + -^^asin(a;j-xOe~"l^^-^'l , (122) 
i=i 

where Xi{t) (E (— oo, +oo) is the opinion of the «th individual at time t and the w^'s are the natural or intrinsic 
opinion changing rates. Note that the interactions are assumed to be all-to-all, though the model can be 
directly generalized to any other topology. Moreover, the w's are drawn from a distribution g{ijj) centered 



at loq with the same properties as in the KM case (Sect. III). In other words, in |255l | the authors approach 
the problem of opinion formation from a radically different point of view in which individuals do not only 
change their opinions, but also the rate at which these changes take place. It is plausible then to assume 
a dynamics described by oscillators coupled together. As opinion changing rates depend on the actual 
interaction between the members of the population, the dynamical model fits quite naively a Kuramoto-like 
model with the additional constraint that individuals having too far opinions will not likely interact. On 
the other hand, as opposite to the KM, opinions are not periodic anymore, so that a new order parameter 
is introduced as 



R{t) 



1 ^ 

\ j^Y.^ij{t)-mf , (123) 

\ 1 



with X{t) being the average of ij{t) over all individuals. From Eq. (|123p . it follows that i? = 1 if complete 
synchronization is achieved and i? < 1 when only partial synchronization occurs. Note that synchronization 
in this context means that the population has reached a unique state of opinion, i.e., it is uniquely polarized 
and that further changes take places at the same rate. Moreover, when complete synchronization is not 
achieved, different opinions are present in the system and the population can be regarded as multipolar. 
The issue is then to investigate under what conditions different emergent behaviors are observed. 

Numerical simulations of the model show that there is a phase transition from incoherence to synchrony 
at a well defined critical coupling CTc. It is argued that when cr < CTc, the society can be thought of as 
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being formed by isolated, non-interacting cultures or groups of opinions, since mixture or agreement is not 
achievable. On the contrary, when <j 3> <Jc, the system fully synchronizes, giving rise in a social context 
to a polarized or globalized society where social and cultural differences are constrained into a single way 
of thinking, notwithstanding the different tendencies to changes of the individuals. Finally, the authors 
reported that bipolarity is only possible if cr ~ Uc, alth ough in this case the model shows a rich behavior 
depending on the way initial opinions are assigned j255| . It would be extremely useful to investigate in the 
future what is the influence of the underlying topology and if the overall picture described above still remains 
valid. Furthermore, as changes in opinions in a population also implies reshaping of the social structure, 
the question of how rewiring mechanisms that take into account the actual opinion states of individuals is 
worth studying in the future. 

5.3.2. Finance 

When reading the economic news, it is not difficult to identify the existence of economic cycles in which 
Gross Domestic Products (GDP's), economic sectors, or stock options raise and fall. Most of the time this 
does not happen for isolated countries, sectors or options but it occurs in quite a synchronized way, although 
some delays are noticeable. 

Within the framework of the current review, we are focusing on synchronization in complex networks, and 
this is what we can identify in many economical sectors: there exists a complicated pattern of interactions 
among companies or countries and the dynamics of each one is quite complex. But, in contrast to many 
networks with a physical background, here we neither know in detail the node dynamics nor its connectivity 
pattern. In this situation it is useful to look at the problem from a different angle. By analyzing some 
macroscopic outcomes, we get some insight into the agents' interactions. 

In the economic literature, synchronization is measured by a correlation coefficient (see, for instance, 

[256| ). based on the idea that correlated (synchronized) business cycles should generate correlated returns. 
The point is to identify what types of interactions lie behind market comovements. Synchronization is 
the result from two different effects. On the one hand, there are different types of common disturbances 
(world interest rates, oil price, or political uncertainty). On the other hand, there exist strong interactions 
between the agents (financial relationships, sector dependencies, co-participation in director boards, etc.). 
It is precisely, these interactions that play a crucial role in the synchronized behavior along economic cycles 
of tightly connected agents and the analysis of the correlations can help in shedding light on the strength 
of the different factors. 

The application of networks co ncepts, mainly that of trees, to economical systems dates back to the 
pioneering work by iMantegna 257 1 , who found a hierarchical arrangement of stocks through the study of 



the correlation returns. 

Recently, the authors in |258| have taken a similar approach to analyze the dynamics of markets. They 
look at the daily closure prices for a total of = 477 stocks traded by the New York Stock Exchange over 
a period of 20 years, from Jan 02, 1980 to Dec 31, 1999. The data is smoothed by looking at time windows 
of given width. As is usually done in the analysis of financial data, the measured quantity is the logarithmic 
return of the stocks, defined as 

n{t)=\nP,{t)-\nP,{t-l), (124) 

where Pi{t) is the closure price of stock i at time t. To quantify the degree of synchronization of the data, 
they use the equal time correlation between assets 

^ {n{t)rM-{r^{t)){r,{t)) ^^^5) 

^[{rm-inim [(r|(t))^(r,(t))2] 

where (. . .) stands for a time average over the consecutive trading days. 

From these correlations the asset tree is constructed. The distance between assets is defined as 



d,,{t) = ^2{1- p,,{t)). (126) 
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The minimum spanning tree is a simply connected graph with iV — 1 edges, such that the sum of aU the 
distances between connected nodes in the graph ^ dij (t) is minimum. 

This procedure generates a time sequence of asset trees that can be interpreted as a sequence of evolu- 
tionary steps of a single dynamics asset tree. For instance, one can identify the leading asset, that, in most 
instances, corresponds to General Electric. 

Such a reduction of the whole set of data retains most of the salient features of the stock market. It is a 
remarkable fact that during crisis periods the markets are very strongly correlated. In terms of the tree its 
average length is reduced and the tree is very tightly packed. By reducing the time window, the location of 
the smallest tree converges to the Black Monday (October 19, 1987). 

It is clear that a hidden pattern of interactions between the assets is responsible for such a synchronized 
behavior and that during crashes the interactions are strengthened. The message here is that the degree 
of synchronization, quantitatively described in terms of the asset correlation, is an indirect measure of the 
existence of strongly connected agents in financial markets. 

The goal of measuring correlations in time series of financial data is to identify synchronized behavior of 
stocks. Stocks are synchronized if they are strongly connected by means of some of the interactions we listed 
above (sharing directors, capital flow, sector dependency, etc.). But clustering thes e stocks according to 
their correlations is usually a hard task and not very accurate. In a recent paper [259j . it has been proposed 
a new method based on synchronization of chaotic maps to get a more precise clustering of the data. They 
look at correlations between stocks in the usual way (c^ Pij{t) >t, is the time average of the correlation 
matrix) but they use these values to construct, through a nonlinear function, a new matrix Jij that is used 
as the interaction matrix between units. The units evolve according to chaotic map dynamics 

Mr + 1) = ^ J^^fi^ji^) (127) 

and /(x) = 1 — 2x^ is the logistic map. The next step is to convert, after some equilibration time, the 
continuous variable Xi into a spin-like one and compute the mutual information between stocks, lij. Then 
stocks are clustered by using lij as the similarity index, and the most stable partition is that with the highest 
cluster entropy. The dendrogram obtained in this way shows, for a particular set of data, a clear partition 
between different classes of stocks. However, the analog procedure applied to the original correlations, Qj, 
shows a "chaining effect" that tends to yield elongated clusters. In this case, synchronization of chaotic 
dynamics turns out to be a powerful tool for the analysis of financial data subjected to similar business 
cycles. 



5.3.3. World Trade Web 

The World Trade Web (WTW) is another example of economic system that has been widely analyzed 
from the network perspective. Different sources provide data about the trade between countries. Networks 
are constructed such that countries are the nodes and trade corresponds to the interaction strengths between 
the nodes . Th ere are several studies that have focused on the static complex natu re of the links between 
countries [260l | and also the evolution of the statistical properties of the network 261 1. But, from a dynamical 
point of view, the trade volume between countries is related to the internal state of the nodes, measured, in 
this case, by the Gross Domestic Product (GDP). 

Hence a clear relationship bet ween node dynamics (which can have a cyclical component) and trade 
strength appears. In particular, in [262| the interplay between the topology of the WTW and the dynamics 
of the GDP's was analyzed. 

In none of the previous studies any reference is made to the precise dynamics of the countries economies. 
As we have stated before, we find what are usually called economic cycles. Economies rise and fall, although 
they mainly raise, but without a constant rate, following rather unpredictable evolutions in time. Due to 
globalization effects, all economies are strongly correlated and they will tend to follow a common trend. In 
our framework, we can say that following a similar time evolution, economies are synch roni zed. This cycle 
synchronization of economies is a topic of current interest in the economic literature; in 263l | a comparative 
analysis between developing and industrial countries is performed, finding that the correlations within the 
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Figure 37: Twenty-two developed countries' economic c ycles synchronization phenomena. The positive real GDP correlation 
means synchronous economic cycles with the US. From [264(1 . 

first group are positive but smaller than within the second group. Another example is found in j264l | where 
correlations between countries are measured in a similar way as the financial series reported above. By 
taking into account that the US is the largest economy, and also the biggest node in the WTW, they look 
at the degree of synchronization between 22 developed countries and the USA (see Fig. 157)) . 

Another clear effect is that particular economies are tightly connected because of economical agreements 
or dependence on particular sectors. In the language of communities this stronger relation can be understood 
as the existence of communities in the overall structure of the world economy. 

Along these lines we have observed that there is a tight relationship between synchronization of economic 
cycles in terms of the GDP and the topological structure of the WTW. The important question raised here 
is if the structure of the network that can be constructed from the empirical data on the cycles correlations, 
as is done with the finance data, can be mapped into the WTW also empirically constructed from the world 
trade transactions. 



6. Perspectives 

As we have seen, the MSF provides a powerful framework to study the relation between the network 
architecture and various types of synchronizability. However, the analysis is mainly limited to the linear 
stability of the complete synchronization states. In most realistic systems where synchronization is relevant, 
complete synchronization of fully identical oscillators is too ideal, and very strong degrees of synchronization 
could relate to pathological activities, such as epileptic seizure in neural systems or social catastrophes. Most 
likely, various levels of synchronization are desirable to enable the system to have flexibility and robustness 
for the emergency of coordination at different scales. 

In this sense, the investigation of synchronization in complex networks is still emerging. We would like 
to list a few future research directions which we think are the main challenges in this field that need to 
be acomplished to achieve a complete comprehension of the structure-synchronization relations in complex 
networks. 



• Spectral properties and synchronization processes 
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In the MSF formulation, we mainly considered the minimal and maximal eigenvalues that are asso- 
ciated to the stability of the synchronization states. The detailed spectral properties, including the 
eigenvectors, of the Laplacian/Adjaceny matrix are important when we are interested in the processes 
leading to synchronization, or interested in the dynamical patterns emerging after a perturbation oc- 
curs. So far studies on deta iled spectral properties are still mainly restricted to random networks 
I23I. [2651 [266L [HI |267| . Ill2 !| and only a few deal on the relation be tween syn chronization patterns 



and spectral properties, out of the complete synchronization, see e.g.. 



• Directed networks and synchronization 

There is growing interest in characterizing directed and weight ed networ ks 269l I60I [l37[[l6ll . [l6l[270t . 



The analysis of the spectra of directed and weighed networks 271 , 272| , which in general are complex, 
and the st udy of the impact of the directed characterizations of the networks on synchronization 
process [9^ llSSL 13], will be key to understanding dynamical organization in more realistic complex 
systems. 

• Co-evolution of structure and synchronization 

Most of the work has considered the impact of network architectures on the synchronization dynamics. 
In many realistic systems, the feedba ck o f dynamics can reshape the network structures, e.g., in neural 
systems through synaptic plasticity |273l |. As a result, adaptation, co-evolution and self-organization 
occur crossing a broad range of scales (see 274l | for a rec ent r e view on adapt ive networks). Adaptation 
due to synchronization has received increasing interest [275I . [27! [nil . [277t . and this line of research 
will lead to an integrated understanding of the structure and dynamics of many complex network 
systems. 



7. Conclusions 

Through the current review we have outlined the state of the art towards a theory of synchronization 
in complex networks. We emphasize the word theory, because, up to now, physicists have made an effort 
of characterization that certainly deepened our understanding of the complex connectivity of natural and 
manmade networks, however, we cannot yet state that we have a theory of complex networks. The topological 
characterization may not be useful to make actual predictions which can be contrasted with experiments. 
To this specific end, the complex network substrate must be enriched and entangled to the functioning of 
the system, i.e., to the dynamics run on top of it. 

The phenomenon of synchronization is one of the paradigmatic observations in different dynamical sys- 
tems. It is at the heart of some biological processes, and according to the wide variety of applications 
presented here, it is a plausible abstraction for many other processes in different contexts. The natural 
approximation to synchronization, from the simplicity of the Kuramoto model, has been explored, and even 
in this case an intricate set of questions concerning the uniformity of the equations proposed, or the nature 
of the critical behavior at the onset of synchronization still have to be definitely settled. The main results, 
however, have helped to understand the nature of the relationship between the topology of interactions and 
the synchronization of phase oscillators. We foresee the importance of these results in the basis of a theory of 
neural dynamics of the brain. Although neural dynamics is far more complex than the phase representation 
reviewed, main features that describe the path towards synchrony in complex networks have been already 
stated and can thus be considered as a good starting point. However, regarding the subjects revised here, 
one easily realizes that the work is far from complete and many questions are in the air, e.g. the evolution 
of synchronization in evolving topologies, the effect on the phenomenon of several kind of disorder, time 
delays, the presence of noise, etc. 

The other main contribution to the problem comes from the MSF formalism. The elegant structure of 
the formalism, designed for linear systems or nonlinear systems close to the synchronization state, allows 
us to make theoretical predictions independent of the specificities of the dynamics. This general framework 
has been deeply studied recently, and provides one of the few mechanisms that allow to make predictions 
about the evolution of synchronized systems as a function of its specific topology. 
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We think that the exploration of new mathematical objects, able to merge the information provided by 
specific dynamics (as for example the Kuramoto model) along to the whole process towards synchronization, 
together with the fine description of the dynamics near the synchronization manifold, should be the focus 
of intense research if we aim to provide a general theory of synchronization processes in complex networks. 
On the other hand, the myriad of applications that can bo cast into mathematical models equivalent to 
those presented along the review, indicates an explosion of activity in different disciplines that will use the 
conceptual framework of synchronization processes in networked systems as a fundamental playground for 
understanding its dynamical behavior. 
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